{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import root,fsolve\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-3319.1198119]\n",
      "[14643.89392131]\n",
      "[-9997.71250365]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([1209.33029053,  587.02356327,  369.35437908,  745.04084551,\n",
       "       1382.24290817, 1587.06989692, 2155.26693601, 1676.77805761,\n",
       "       2740.87800159,  817.80899893, 3512.46974556,  109.55037279])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(1)\n",
    "    #nois_dist = np.random.\n",
    "    #axis = np.random.uniform(low,high,(num,3))\n",
    "print(np.random.uniform(-20000,20000,size=(1,)))\n",
    "print(np.random.uniform(6000,18000,size =(1,)))\n",
    "print(np.random.uniform(-10000,10000,size =(1,)))\n",
    "np.random.uniform(0,4000,size =(12,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 2., 3.],\n",
       "       [1., 2., 3.],\n",
       "       [1., 2., 3.]])"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array((1,2,3)).reshape(1,3) -np.zeros((3,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.00315281, -0.00253832, -0.00627945, -0.00219859, -0.00139763,\n",
       "        0.00682319, -0.00424491, -0.0020083 , -0.00300981, -0.00607778,\n",
       "        0.0004164 , -0.00062512])"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.normal(0,dX,(num,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "np.random.normal(0,dL,(num,));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "8*np.arange(0,7,1) +4;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4 [8.60313976 4.06179954 2.82596356] 9.924630940063436 9.513791529699938\n",
      "8 [4.79135136 7.80404596 2.02323567] 9.378361463371865 9.157520470077603\n",
      "12 [ 4.15107392 12.92085584  1.70861166] 13.678424037668927 13.571290666402302\n",
      "16 [ 3.0882686  12.18797322  1.14256325] 12.624957222267682 12.573149727790131\n"
     ]
    }
   ],
   "source": [
    "#print(myint) 正侧\n",
    "iters =20000\n",
    "angle = 20\n",
    "y_axis = 0\n",
    "z_axis = 10000\n",
    "x_axis = 0\n",
    "\n",
    "numGrp = 4*np.arange(1,5,1) #+4#2*np.arange(2,24,2)\n",
    "#numGrp=[4]\n",
    "result = []\n",
    "for num in numGrp:\n",
    "    calu = np.zeros(3)\n",
    "    for i in range(iters):\n",
    "        myint = np.random.randint(0,20000)\n",
    "        np.random.seed(myint)\n",
    "\n",
    "        Paxis = np.zeros((num,3))\n",
    "        Daxis = np.zeros((num,3))\n",
    "\n",
    "        xyz0 = np.array((x_axis,y_axis,z_axis)).squeeze()\n",
    "\n",
    "        dX = 0.03*100\n",
    "        dY = 0.03*100\n",
    "        dZ = 0.05*100\n",
    "        dL = 0.02*100  #0.02\n",
    "        # print(x_axis)\n",
    "        # print(y_axis)\n",
    "        # print(z_axis)\n",
    "        Daxis[:,0] = np.random.normal(0,dX,(num,))\n",
    "        Daxis[:,1] = np.random.normal(0,dY,(num,))\n",
    "        Daxis[:,2] = np.random.normal(0,dZ,(num,))\n",
    "        DL = np.random.normal(0,dL,(num,))\n",
    "        #print(DL)\n",
    "\n",
    "        dx0 = z_axis*np.tan(angle*np.pi/180)\n",
    "        dx1 = dx0 + 1000 #5000 \n",
    "        #print(dx0, dx1)\n",
    "\n",
    "        #dz0 = z_axis*np.tan(20*np.pi/180)\n",
    "        #dz1 =dz0 + 2500\n",
    "        \n",
    "        Paxis[:,0] = np.random.uniform(dx0, dx1,size =(num,))\n",
    "        #print(Paxis)\n",
    "        Paxis[:,1] = 0 #np.random.uniform(-1500, 1500,size =(num,)) #0\n",
    "        Paxis[:,2] = 0 #np.random.uniform(0, 0,size =(num,))  \n",
    "        #Paxis[:num//2+1,0] = -Paxis[:num//2+1,0]\n",
    "        Paxis[:num//2,0] = -Paxis[:num//2,0]\n",
    "        \n",
    "        distL0 = np.sqrt( (x_axis-Paxis[:,0])**2 + (y_axis-Paxis[:,1])**2+ (z_axis-Paxis[:,2])**2 )\n",
    "        #print(distL0)\n",
    "        distL1 = distL0 + DL\n",
    "        distL1 = distL1.reshape(num,1)\n",
    "        #print(distL1)\n",
    "        ##########################\n",
    "        #print(Paxis)\n",
    "        Paxis += Daxis\n",
    "        #print(Paxis)\n",
    "        results = []\n",
    "\n",
    "\n",
    "        def f5(x):\n",
    "            m = np.array((x[0],x[1],x[2]),dtype = np.float64).reshape(1,3)\n",
    "            minus = (m - Paxis)\n",
    "        #    print('###################')\n",
    "        #    print(Paxis)\n",
    "        #    print('#####################')\n",
    "            k = np.concatenate((minus**2, -distL1**2),axis = 1)\n",
    "            l = np.sum(k,keepdims = True,axis =1).reshape(-1,1)\n",
    "            return np.sum(l * minus,axis = 0)\n",
    "\n",
    "        #Sxyz_tmp = root(f5,[x_axis,y_axis,z_axis])\n",
    "\n",
    "        Sxyz = fsolve(f5,[x_axis,y_axis,z_axis])\n",
    "\n",
    "        #D1xyz = Sxyz_tmp.x -xyz0\n",
    "        Dxyz = (Sxyz - xyz0)**2\n",
    "        calu += Dxyz\n",
    "        #print(xyz0)\n",
    "        #print(Sxyz)\n",
    "    calu /= iters\n",
    "    #print(np.concatenate((np.array((num,)),np.sqrt(calu), np.sqrt(np.sum(calu)),np.sqrt(np.sum(calu[0:2])))))\n",
    "    result.append(num)\n",
    "    result += np.sqrt(calu).tolist()\n",
    "    result.append(np.sqrt(np.sum(calu)))\n",
    "    result.append(np.sqrt(np.sum(calu[0:2])))\n",
    "    print(num, np.sqrt(calu), np.sqrt(np.sum(calu)),np.sqrt(np.sum(calu[0:2])))\n",
    "        #print(D1xyz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ToExcel(result,file_path):\n",
    "    mode = 'a' if os.path.exists(file_path) else 'w'\n",
    "    writer = pd.ExcelWriter(file_path,mode=mode)\n",
    "    result.to_excel(writer, index=False,encoding='utf-8',sheet_name='Sheet')\n",
    "    writer.save()\n",
    "result = pd.DataFrame(np.array(result).reshape(-1,6),columns=['NUM','σ_X','σ_Y','σ_Z','σ_R','σ_XY'],)\n",
    "file_path = 'C://Users//gao//Desktop//1.xlsx'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "ToExcel(result,file_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4 [130.76304654 127.34952949  10.80689804] 182.8487518277908 182.5291127454484\n",
      "12 [30.3943958  29.76667092  3.19269208] 42.66224650292494 42.542613859090885\n",
      "20 [24.3866226  23.56347426  2.48582245] 34.001823394107774 33.910834269495\n",
      "28 [18.41565761 17.6068813   1.57394607] 25.526770664166186 25.4782007668851\n",
      "36 [16.85366794 16.41990214  1.77831078] 23.597069705150258 23.529966201296595\n",
      "44 [16.58783984 15.76916437  1.73032023] 22.952494059209265 22.887179280992925\n",
      "52 [12.99272254 12.93636705  1.25737055] 18.37774229803184 18.33467837948326\n"
     ]
    }
   ],
   "source": [
    "######################前左+后右##################\n",
    "iters =20000\n",
    "angle = 40\n",
    "y_axis = 0\n",
    "z_axis = 10000\n",
    "x_axis = 0\n",
    "\n",
    "numGrp = 8*np.arange(0,7,1) +4#2*np.arange(2,24,2)\n",
    "#numGrp=[4]\n",
    "result = []\n",
    "for num in numGrp:\n",
    "    calu = np.zeros(3)\n",
    "    for i in range(iters):\n",
    "        myint = np.random.randint(0,20000)\n",
    "        np.random.seed(myint)\n",
    "\n",
    "        Paxis = np.zeros((num,3))\n",
    "        Daxis = np.zeros((num,3))\n",
    "\n",
    "        xyz0 = np.array((x_axis,y_axis,z_axis)).squeeze()\n",
    "\n",
    "        dX = 0.05*100\n",
    "        dY = 0.05*100\n",
    "        dZ = 0.05*100\n",
    "        dL = 0.03*100  #0.02\n",
    "        # print(x_axis)\n",
    "        # print(y_axis)\n",
    "        # print(z_axis)\n",
    "        Daxis[:,0] = np.random.normal(0,dX,(num,))\n",
    "        Daxis[:,1] = np.random.normal(0,dY,(num,))\n",
    "        Daxis[:,2] = np.random.normal(0,dZ,(num,))\n",
    "        DL = np.random.normal(0,dL,(num,))\n",
    "        #print(DL)\n",
    "\n",
    "        dx0 = z_axis*np.tan(angle*np.pi/180)\n",
    "        dx1 = dx0 + 2500 #2500\n",
    "        #print(dx0, dx1)\n",
    "\n",
    "        dz0 = z_axis*np.tan(40*np.pi/180)\n",
    "        dz1 =dz0 +3000\n",
    "        \n",
    "        Paxis[:,0] = np.random.uniform(dx0, dx1,size =(num,))\n",
    "        #print(Paxis)\n",
    "        Paxis[:num,1] = -np.random.uniform(dz0, dz1,size =(num,))#np.random.uniform(0, 3000,size =(num,)) \n",
    "        #0 #np.random.uniform(dz0, dz1,size =(num//2,))\n",
    "        Paxis[:,2] = 0 #np.random.uniform(0, 0,size =(num,))  \n",
    "        #Paxis[:num//2+1,0] = -Paxis[:num//2+1,0]\n",
    "        Paxis[:num//2,0] = -Paxis[:num//2,0]\n",
    "        Paxis[:num//2,1] = -Paxis[:num//2,1]\n",
    "        \n",
    "        distL0 = np.sqrt( (x_axis-Paxis[:,0])**2 + (y_axis-Paxis[:,1])**2+ (z_axis-Paxis[:,2])**2 )\n",
    "        #print(distL0)\n",
    "        distL1 = distL0 + DL\n",
    "        distL1 = distL1.reshape(num,1)\n",
    "        #print(distL1)\n",
    "        ##########################\n",
    "        #print(Paxis)\n",
    "        Paxis += Daxis\n",
    "        #print(Paxis)\n",
    "        results = []\n",
    "\n",
    "\n",
    "        def f5(x):\n",
    "            m = np.array((x[0],x[1],x[2]),dtype = np.float64).reshape(1,3)\n",
    "            minus = (m - Paxis)\n",
    "        #    print('###################')\n",
    "        #    print(Paxis)\n",
    "        #    print('#####################')\n",
    "            k = np.concatenate((minus**2, -distL1**2),axis = 1)\n",
    "            l = np.sum(k,keepdims = True,axis =1).reshape(-1,1)\n",
    "            return np.sum(l * minus,axis = 0)\n",
    "\n",
    "        #Sxyz_tmp = root(f5,[x_axis,y_axis,z_axis])\n",
    "\n",
    "        Sxyz = fsolve(f5,[x_axis,y_axis,z_axis])\n",
    "\n",
    "        #D1xyz = Sxyz_tmp.x -xyz0\n",
    "        Dxyz = (Sxyz - xyz0)**2\n",
    "        calu += Dxyz\n",
    "        #print(xyz0)\n",
    "        #print(Sxyz)\n",
    "    calu /= iters\n",
    "    #print(np.concatenate((np.array((num,)),np.sqrt(calu), np.sqrt(np.sum(calu)),np.sqrt(np.sum(calu[0:2])))))\n",
    "    result.append(num)\n",
    "    result += np.sqrt(calu).tolist()\n",
    "    result.append(np.sqrt(np.sum(calu)))\n",
    "    result.append(np.sqrt(np.sum(calu[0:2])))\n",
    "    print(num, np.sqrt(calu), np.sqrt(np.sum(calu)),np.sqrt(np.sum(calu[0:2])))\n",
    "        #print(D1xyz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = pd.DataFrame(np.array(result).reshape(-1,6),columns=['NUM','σ_X','σ_Y','σ_Z','σ_R','σ_XY'],)\n",
    "file_path = 'C://Users//gao//Desktop//1.xlsx'\n",
    "writer = pd.ExcelWriter(file_path)\n",
    "result.to_excel(writer, index=False,encoding='utf-8',sheet_name='Sheet')\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4 [ 6.44461736 68.07174584  6.76657068] 68.71013136957556 68.37613380451013\n",
      "12 [ 2.32431934 13.81771189  1.79927873] 14.126890181832973 14.01183864681278\n",
      "20 [1.5935153  8.26065664 1.52066066] 8.549277626613279 8.412950676392862\n",
      "28 [2.01600875 2.42179835 1.14460135] 3.3525379604681493 3.151094814125687\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-20-e78b7d269ca4>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     30\u001b[0m         \u001b[0mDaxis\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnormal\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdY\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnum\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     31\u001b[0m         \u001b[0mDaxis\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnormal\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdZ\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnum\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 32\u001b[1;33m         \u001b[0mDL\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnormal\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdL\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnum\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     33\u001b[0m         \u001b[1;31m#print(DL)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     34\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "###########前后#######\n",
    "iters =20000\n",
    "angle = 40\n",
    "y_axis = 0\n",
    "z_axis = 10000\n",
    "x_axis = 0\n",
    "\n",
    "numGrp = 8*np.arange(0,7,1) +4#2*np.arange(2,24,2)\n",
    "#numGrp=[4]\n",
    "result = []\n",
    "for num in numGrp:\n",
    "    calu = np.zeros(3)\n",
    "    for i in range(iters):\n",
    "        myint = np.random.randint(0,20000)\n",
    "        np.random.seed(myint)\n",
    "\n",
    "        Paxis = np.zeros((num,3))\n",
    "        Daxis = np.zeros((num,3))\n",
    "\n",
    "        xyz0 = np.array((x_axis,y_axis,z_axis)).squeeze()\n",
    "\n",
    "        dX = 0.03*100\n",
    "        dY = 0.03*100\n",
    "        dZ = 0.05*100\n",
    "        dL = 0.02*100  #0.02\n",
    "        # print(x_axis)\n",
    "        # print(y_axis)\n",
    "        # print(z_axis)\n",
    "        Daxis[:,0] = np.random.normal(0,dX,(num,))\n",
    "        Daxis[:,1] = np.random.normal(0,dY,(num,))\n",
    "        Daxis[:,2] = np.random.normal(0,dZ,(num,))\n",
    "        DL = np.random.normal(0,dL,(num,))\n",
    "        #print(DL)\n",
    "\n",
    "        dx0 = z_axis*np.tan(angle*np.pi/180)\n",
    "        dx1 = dx0 + 1000 #2500\n",
    "        #print(dx0, dx1)\n",
    "\n",
    "        #dz0 = z_axis*np.tan(40*np.pi/180)\n",
    "        #dz1 =dz0 +3000\n",
    "        \n",
    "        Paxis[:,0] = np.random.uniform(dx0, dx1,size =(num,))\n",
    "        #print(Paxis)\n",
    "        Paxis[:,1] = np.random.uniform(-1000, 1000,size =(num,))#np.random.uniform(0, 3000,size =(num,)) \n",
    "        #0 #np.random.uniform(dz0, dz1,size =(num//2,))\n",
    "        Paxis[:,2] = 0 #np.random.uniform(0, 0,size =(num,))  \n",
    "        #Paxis[:num//2+1,0] = -Paxis[:num//2+1,0]\n",
    "        Paxis[:num//2,0] = -Paxis[:num//2,0]\n",
    "        \n",
    "        distL0 = np.sqrt( (x_axis-Paxis[:,0])**2 + (y_axis-Paxis[:,1])**2+ (z_axis-Paxis[:,2])**2 )\n",
    "        #print(distL0)\n",
    "        distL1 = distL0 + DL\n",
    "        distL1 = distL1.reshape(num,1)\n",
    "        #print(distL1)\n",
    "        ##########################\n",
    "        #print(Paxis)\n",
    "        Paxis += Daxis\n",
    "        #print(Paxis)\n",
    "        results = []\n",
    "\n",
    "\n",
    "        def f5(x):\n",
    "            m = np.array((x[0],x[1],x[2]),dtype = np.float64).reshape(1,3)\n",
    "            minus = (m - Paxis)\n",
    "        #    print('###################')\n",
    "        #    print(Paxis)\n",
    "        #    print('#####################')\n",
    "            k = np.concatenate((minus**2, -distL1**2),axis = 1)\n",
    "            l = np.sum(k,keepdims = True,axis =1).reshape(-1,1)\n",
    "            return np.sum(l * minus,axis = 0)\n",
    "\n",
    "        #Sxyz_tmp = root(f5,[x_axis,y_axis,z_axis])\n",
    "\n",
    "        Sxyz = fsolve(f5,[x_axis,y_axis,z_axis])\n",
    "\n",
    "        #D1xyz = Sxyz_tmp.x -xyz0\n",
    "        Dxyz = (Sxyz - xyz0)**2\n",
    "        calu += Dxyz\n",
    "        #print(xyz0)\n",
    "        #print(Sxyz)\n",
    "    calu /= iters\n",
    "    #print(np.concatenate((np.array((num,)),np.sqrt(calu), np.sqrt(np.sum(calu)),np.sqrt(np.sum(calu[0:2])))))\n",
    "    result.append(num)\n",
    "    result += np.sqrt(calu).tolist()\n",
    "    result.append(np.sqrt(np.sum(calu)))\n",
    "    result.append(np.sqrt(np.sum(calu[0:2])))\n",
    "    print(num, np.sqrt(calu), np.sqrt(np.sum(calu)),np.sqrt(np.sum(calu[0:2])))\n",
    "        #print(D1xyz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = pd.DataFrame(np.array(result).reshape(-1,6),columns=['NUM','σ_X','σ_Y','σ_Z','σ_R','σ_XY'],)\n",
    "file_path = 'C://Users//gao//Desktop//1.xlsx'\n",
    "writer = pd.ExcelWriter(file_path)\n",
    "result.to_excel(writer, index=False,encoding='utf-8',sheet_name='Sheet')\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4 [114.13341029 255.88362549 107.86266575] 300.2286125642999 280.1836275373938\n",
      "12 [32.89677429 71.76424535 32.12734656] 85.23186649184632 78.9449470743959\n",
      "20 [30.8008204  67.27915793 30.3384813 ] 79.97248950162435 73.99442971992025\n",
      "28 [21.88367708 47.01725313 21.02236756] 55.95942594405107 51.860557403279216\n",
      "36 [15.22178942 33.15575357 15.1751695 ] 39.5131957366685 36.482966821495275\n",
      "44 [14.44132625 31.2362585  12.90134208] 36.751875826347444 34.413017148562744\n",
      "52 [14.97093311 32.23092335 14.29347813] 38.304892319074796 35.53816622347996\n"
     ]
    }
   ],
   "source": [
    "###########正侧+前斜#######\n",
    "iters =20000\n",
    "angle = 40\n",
    "y_axis = 0\n",
    "z_axis = 10000\n",
    "x_axis = 0\n",
    "\n",
    "numGrp = 8*np.arange(0,7,1) +4#2*np.arange(2,24,2)\n",
    "#numGrp=[4]\n",
    "result = []\n",
    "for num in numGrp:\n",
    "    calu = np.zeros(3)\n",
    "    for i in range(iters):\n",
    "        myint = np.random.randint(0,20000)\n",
    "        np.random.seed(myint)\n",
    "\n",
    "        Paxis = np.zeros((num,3))\n",
    "        Daxis = np.zeros((num,3))\n",
    "\n",
    "        xyz0 = np.array((x_axis,y_axis,z_axis)).squeeze()\n",
    "\n",
    "        dX = 0.05*100\n",
    "        dY = 0.05*100\n",
    "        dZ = 0.05*100\n",
    "        dL = 0.03*100  #0.02\n",
    "        # print(x_axis)\n",
    "        # print(y_axis)\n",
    "        # print(z_axis)\n",
    "        Daxis[:,0] = np.random.normal(0,dX,(num,))\n",
    "        Daxis[:,1] = np.random.normal(0,dY,(num,))\n",
    "        Daxis[:,2] = np.random.normal(0,dZ,(num,))\n",
    "        DL = np.random.normal(0,dL,(num,))\n",
    "        #print(DL)\n",
    "\n",
    "        dx0 = z_axis*np.tan(angle*np.pi/180)\n",
    "        dx1 = dx0 + 2500 #2500\n",
    "        #print(dx0, dx1)\n",
    "\n",
    "        dz0 = z_axis*np.tan(40*np.pi/180)\n",
    "        dz1 =dz0 +1000\n",
    "        \n",
    "        Paxis[:,0] = np.random.uniform(dx0, dx1,size =(num,))\n",
    "        #print(Paxis)\n",
    "        Paxis[:num//2,1] = np.random.uniform(dz0, dz1,size =(num//2,))#np.random.uniform(0, 3000,size =(num,)) \n",
    "        #0 #np.random.uniform(dz0, dz1,size =(num//2,))\n",
    "        Paxis[:,2] = 0 #np.random.uniform(0, 0,size =(num,))  \n",
    "        #Paxis[:num//2+1,0] = -Paxis[:num//2+1,0]\n",
    "        Paxis[:num//2,0] = -Paxis[:num//2,0]\n",
    "        \n",
    "        distL0 = np.sqrt( (x_axis-Paxis[:,0])**2 + (y_axis-Paxis[:,1])**2+ (z_axis-Paxis[:,2])**2 )\n",
    "        #print(distL0)\n",
    "        distL1 = distL0 + DL\n",
    "        distL1 = distL1.reshape(num,1)\n",
    "        #print(distL1)\n",
    "        ##########################\n",
    "        #print(Paxis)\n",
    "        Paxis += Daxis\n",
    "        #print(Paxis)\n",
    "        results = []\n",
    "\n",
    "\n",
    "        def f5(x):\n",
    "            m = np.array((x[0],x[1],x[2]),dtype = np.float64).reshape(1,3)\n",
    "            minus = (m - Paxis)\n",
    "        #    print('###################')\n",
    "        #    print(Paxis)\n",
    "        #    print('#####################')\n",
    "            k = np.concatenate((minus**2, -distL1**2),axis = 1)\n",
    "            l = np.sum(k,keepdims = True,axis =1).reshape(-1,1)\n",
    "            return np.sum(l * minus,axis = 0)\n",
    "\n",
    "        #Sxyz_tmp = root(f5,[x_axis,y_axis,z_axis])\n",
    "\n",
    "        Sxyz = fsolve(f5,[x_axis,y_axis,z_axis])\n",
    "\n",
    "        #D1xyz = Sxyz_tmp.x -xyz0\n",
    "        Dxyz = (Sxyz - xyz0)**2\n",
    "        calu += Dxyz\n",
    "        #print(xyz0)\n",
    "        #print(Sxyz)\n",
    "    calu /= iters\n",
    "    #print(np.concatenate((np.array((num,)),np.sqrt(calu), np.sqrt(np.sum(calu)),np.sqrt(np.sum(calu[0:2])))))\n",
    "    result.append(num)\n",
    "    result += np.sqrt(calu).tolist()\n",
    "    result.append(np.sqrt(np.sum(calu)))\n",
    "    result.append(np.sqrt(np.sum(calu[0:2])))\n",
    "    print(num, np.sqrt(calu), np.sqrt(np.sum(calu)),np.sqrt(np.sum(calu[0:2])))\n",
    "        #print(D1xyz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = pd.DataFrame(np.array(result).reshape(-1,6),columns=['NUM','σ_X','σ_Y','σ_Z','σ_R','σ_XY'],)\n",
    "file_path = 'C://Users//gao//Desktop//1.xlsx'\n",
    "writer = pd.ExcelWriter(file_path)\n",
    "result.to_excel(writer, index=False,encoding='utf-8',sheet_name='Sheet')\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3, 6)\n",
      "[[-5.05207704e+08]\n",
      " [-4.81862312e+08]\n",
      " [-4.98217660e+08]]\n",
      "(3, 6)\n",
      "[[-5.04770168e+08]\n",
      " [-4.81435196e+08]\n",
      " [-4.97783200e+08]]\n"
     ]
    }
   ],
   "source": [
    "print(f6(Sxyz))\n",
    "print(f6(xyz0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndentationError",
     "evalue": "unexpected indent (<ipython-input-119-d7280cd94844>, line 2)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-119-d7280cd94844>\"\u001b[1;36m, line \u001b[1;32m2\u001b[0m\n\u001b[1;33m    results.append(np.sqrt(calu))\u001b[0m\n\u001b[1;37m    ^\u001b[0m\n\u001b[1;31mIndentationError\u001b[0m\u001b[1;31m:\u001b[0m unexpected indent\n"
     ]
    }
   ],
   "source": [
    "   \n",
    "    calu += (Sxyz - xyz0)**2\n",
    "        results.append(np.sqrt(calu))\n",
    "    results = np.array(results).reshape(-1,3)\n",
    "    \n",
    "    range_x1 = np.linspace(1,3*3,17)\n",
    "    plt.figure(k)\n",
    "    plt.title(\"%d\" %(num))\n",
    "    plt.plot(range_x1,results[:,0],'r',range_x1,results[:,1],'b',range_x1,results[:,2],'g')\n",
    "    print(k,results[5,0],results[5,1],results[5,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.26794919243112275"
      ]
     },
     "execution_count": 120,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(3-3**0.5)/(3+3**0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[7182.40479878, 3907.26339744, 2922.96107164],\n",
       "       [7776.16831892, 3907.26339744, 2314.95756424],\n",
       "       [6858.63023588, 3907.26339744, 3236.07067848],\n",
       "       [9434.69548527, 3907.26339744,  668.9940166 ],\n",
       "       [8845.56673495, 3907.26339744, 1248.57987431],\n",
       "       [8666.13046671, 3907.26339744, 1431.72621587],\n",
       "       [8125.72088085, 3907.26339744, 1966.98919628],\n",
       "       [8129.30327613, 3907.26339744, 1960.79295835]])"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[11188.88178083 11796.88528822 10875.77217398 13442.84883587]\n"
     ]
    }
   ],
   "source": [
    "num = 4\n",
    "np.random.seed(102)\n",
    "    #nois_dist = np.random.\n",
    "    #axis = np.random.uniform(low,high,(num,3))\n",
    "axis = np.zeros((num,3))\n",
    "y_axis = np.random.uniform(-20000,20000,size=(1,))\n",
    "z_axis = np.random.uniform(6000,18000,size =(1,))\n",
    "x_axis = np.random.uniform(-10000,10000,size =(1,))\n",
    "   # print(x_axis)\n",
    "   # print(y_axis)\n",
    "   # print(z_axis)\n",
    "axis[:,1] = y_axis\n",
    "axis[:,2] = np.random.uniform(0,4000,size =(num,))\n",
    "b = np.abs(z_axis-axis[:,2])\n",
    "print(b)\n",
    "axis[:,0] = np.random.uniform(6000,8500,size =(num,))# x_axis + np.abs(z_axis-axis[:,2]) +np.random.uniform(-10,10,(num,))\n",
    "\n",
    "#np.random.uniform((z_axis-axis[:,2]) /3**0.5,(z_axis-axis[:,2]) * 3**0.5)\n",
    "#np.random.uniform((x_axis-1,x_axis+1),3)\n",
    "gt = np.array((x_axis,y_axis,z_axis)).squeeze()\n",
    "#slope_distance= np.random.randint(10,20,cal).reshape(cal,1)\n",
    "slope_distance = np.sqrt((x_axis-axis[:,0])**2 + (z_axis-axis[:,2])**2) .reshape(num,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[6780.36242144, 3907.26339744, 2922.96107164],\n",
       "       [6894.82888492, 3907.26339744, 2314.95756424],\n",
       "       [7229.36824768, 3907.26339744, 3236.07067848],\n",
       "       [7225.49559897, 3907.26339744,  668.9940166 ]])"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x20b7bb2e148>,\n",
       " <matplotlib.lines.Line2D at 0x20b7bdf5108>,\n",
       " <matplotlib.lines.Line2D at 0x20b7be01d08>]"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEICAYAAABGaK+TAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZyN5f/H8dc1MxgG2UbWQllTyJCQr7JUikK2/CRr30JShEgLifJVCGUnkT0SkiUUydjKEmXfTdmXMWbm+v1xH8swMWPOuM+ZeT8fj/M4c+5zn3N/jsd4u1znWoy1FhER8T8BbhcgIiK3RgEuIuKnFOAiIn5KAS4i4qcU4CIifkoBLiLipxTgIiJ+SgEuAhhjChtjIo0xE92uRSShFOAijqHAGreLEEkMBbikesaYxsAJYLHbtYgkhgJcUjVjTGbgfeANt2sRSSwFuKR2vYHR1tp9bhciklhBbhcg4hZjTGmgOlDG7VpEboUCXFKzqkABYK8xBiAjEGiMKWGtfdDFukQSxGg5WUmtjDEZgMxXHeqME+gvW2sjXClKJBHUApdUy1p7Djh36bEx5gwQqfAWf6EWuIiIn9IoFBERP6UAFxHxUwpwERE/pQAXEfFTt3UUSo4cOWyBAgVu5yVFRPze2rVr/7bWhl57/LYGeIECBQgPD7+dlxQR8XvGmD3xHVcXioiIn1KAi4j4KQW4iIifUoCLiPgpBbiIiJ9SgIuI+CkFuIiIn1KAi4gkp6NH4fXX4fRpr7+1AlxEJLlYCy+/DEOHwt69Xn97beggIpJcJk+GmTOhf3+47z6vv71a4CIiyeHgQWjfHh5+GN54I1kuoQAXEfE2a6FNG4iMhHHjIDAwWS6jLhQREW8bOxbmzYNBg6BIkWS7jFrgIiLetGcPvPYaVK3qdKEkIwW4iIi3xMZCq1ZOF8qYMRCQvBGrLhQREW/5/HNYvBi++AIKFkz2y6kFLiLiDTt2QJcu8PjjzheYt4ECXEQkqWJi4MUXIU0aGDUKjLktl1UXiohIUg0aBD/9BOPHQ758t+2yCWqBG2M6GWM2G2M2GWMmG2OCjTEFjTGrjTF/GmOmGGPSJnexIiI+Z+tWeOstqFMHmjW7rZe+aYAbY/ICrwJh1tqSQCDQGOgPfGKtLQwcB1olZ6EiIj4nOhqaN4eMGZ0vLuPpOjly5gjtvmvH6QvuLWYVBKQ3xgQBGYBDwGPAdM/z44FnvV6diIgv++gjWLMGhg2DXLmuezomNobnZz7P2A1j2XvS+4tZ3TTArbUHgAHAXpzgPgmsBU5Ya6M9p+0H8sb3emNMW2NMuDEmPCIiwjtVi4i47bff4N13oWFD5xaP95e9z5JdSxhaayj35XRhMStjTFbgGaAgkAcIAZ6M51Qb3+uttSOstWHW2rDQ0NCk1Coi8u/OnIFp0yAqKvmvFRUFL7wA2bI5S8XG44cdP9B7eW+al2pOizItkqWMhHShVAd2WWsjrLUXgZlARSCLp0sFIB9wMFkqFBFJiO7dnZZwmTKwcmXyXqtPH9i4EUaMgBw5rnv64OmDNJ3ZlBKhJRhaK/6A94aEBPheoIIxJoMxxgDVgC3AUuA5zznNgdnJU6KIyE3s3euEabVqTku8UiV45RU4edL71woPh759nS8v69S57uno2GgaT2/MuYvnmNZgGiFpQ7xfg0dC+sBX43xZuQ743fOaEUBX4HVjzF9AdmB0slUpInIjH3zg3I8ZA5s3Q6dOzqiQ4sVhxgxnbRJviIx0uk5y5YJPP433lF5Le7Fi7wo+f/pziocW9851/4219rbdypYta0VEvGrHDmuDgqxt1y7u8TVrrC1d2lqwtk4da/fuTfq1unRx3m/Bgnif/m77d5Z3sW3mtEn6ta4ChNt4MlVT6UXEv/XuDUFBzmSaq4WFOUP8BgyARYugRAkYMsSZ9n4rVq503uull5z1Tq6x7+Q+ms1qRqk7SzHoiUG3do1EUoCLiP/avh0mTHA2Ds6T5/rng4Kc7cw2bYLKleHVV6FiRecLyMQ4e9bp8777bvj44+uevhhzkUbTGxEVE8XUBlNJnyb9LX6gxFGAi4j/eu89CA6Gbt1ufF7Bgs4OOZMmwa5dULas85pz5xJ2ne7d4a+/nJ12MmW6/unF3Vm1fxWjao+iSPbk24HnWgpwEfFPmzc7u7536AA5c978fGOgSRP44w9n5cD+/aFkSVi48MavW7rU6Xrp2NHZZecac7bN4X+r/scrYa/QqGSjW/oot8pYb307mwBhYWE2PDz8tl1PRFKwBg1gwQKnRR3PWOyb+vFHpz97+3b4v/+DgQPh2smGp07BAw9A2rSwYQNkyBDn6d0ndlPmizIUylqIlS1Xki4o3a1/nhswxqy11oZde1wtcBHxPxs3wvTpzt6TtxLe4LSmN26Et9+GKVOgWDFnOdirG7WdO8O+fc7xa8I7KiaKhtMaYq1lWoNpyRbeN6IAFxH/8847cMcd8Prrlw+dvnCaKZumEBkdmfD3CQ6G9993WtfFiztdK9Wrw59/Oq37kSOdXXYefvi6l3ZZ2IU1B9cw9pmxFMpayAsfKvEU4CLiX8LDYfZsZ3RJ1qyXD7+x8A0az2jMfcPu45s/viFR3cMlSsDy5c6eluHhcP/9TrfKffc5X5ReY8aWGQz+dTCvPfQadYvX9canuiUKcBHxL716OYtIdex4+dCmo5sYvX40zxR9huCgYOpOqcvjEx9nS8SWhL9vQIDTJ751K9SuDefPO10n6eJ2jew4toOWc1pSPm95+tfo761PdUsU4CLiP1atgvnz4c03IXPmy4ff/OFNMqfLzOg6o9nw0gYGPTGINQfX8MDwB3htwWuciDyR8GvkyeOsanjihDPc8CqR0ZE0mNaAQBPI1OemkjbQ3Y3IFOAi4j/eftsZKdK+/eVDP+z4gfl/zafnIz3JniE7aQLT8OpDr7K9/XZalWnF4NWDKTykMCPXjiQmNhGzMNOkue7Q69+/zvrD65lQdwJ3Z7nbG58oSRTgIuIfli2DxYudCTghzgp/MbExdP6hMwWyFKB9+fZxTg8NCeWL2l+wtu1aiuUoRtu5bSk3shw/7/35li4/+ffJDA8fTpeKXXi6yNNJ/jjeoAAXEd9nrdP6zp3bmTbv8eVvX/Lbkd/oV63fvw7jK5O7DMtfXM7k+pM5evYolcdWpunMpuw/tT/Bl9/29zbazm1LxfwV+eCxD5L8cbxFAS4ivm/RIlixwlmwKr2zzsjZqLP0WNKDh/I+RMP74t/S7BJjDI1LNmZb+230fKQnM7bMoOhnRflg+Qc3HXZ47uI5GkxrQLrAdEx5bgppAq/vWnGLAlxEfNul1nf+/NCmzeXDA1cN5ODpg/yv5v8w8ewGH5+QtCH0fqw3W9tt5fF7Hqfn0p6UGFrihsMOX53/Kr8f/Z2J9SaSL3M+r3wkb1GAi4hvmz8fVq+Gnj0vD+k7fOYw/X/uT/3i9al0V6VEv2XBrAWZ2Wgmi5otIkOaDNSdUpeaE2teN+xwwsYJjF4/mrcqv8UT9z7hlY/jTVoLRUR8l7VQrhwcOwbbtl0eGfLSty8xZsMYtrbbyr3Z7k3SJaJjoxm+Zji9fuzF6QunaV++Pe9WfZeDpw9SbmQ5yuUpx6IXFhEUEHTzN0sm/7YWinsViYjczOzZsHats4yrJ7w3H93MqPWj6FC+Q5LDGyAoIIgOD3WgccnGvL30bQavHsxXv39FprSZyJg2I5PrT3Y1vG9EXSgi4ptiY51Zl4ULO9PaPd5c9CaZ0mbi7Spve/VyoSGhfP7055eHHe49uZdJ9SaRO1Nur17Hm3zznxURkenT4fffYeJEZ2cdYNHORcz7cx4f1/iY7BmyJ8tlLw07PB55nGzpsyXLNbxFfeAi4ntiYpwFpYyB336DwEBiYmMoO6IsJy+cZGu7rQQHBbtd5W2jPnAR8R+TJzuLSk2dCoGBgDNpZ+ORjUyuPzlVhfeNqAUuIr4lOtpZmztDBli/HgICOHfxHIWHFCZf5nz80uqXBI/7TinUAhcR/zBhgrOB8DffOEu8cmXSztf1v0514X0jGoUiIr4jKgp694awMKhTB3Am7fT7qR91i9XlkbsfcblA36IWuIj4jrFjYfduGDbM+QITePfHd7kQc4F+1fu5W5sPUgtcRHxDZCT06QMVKsATzrT1LRFbGLluJC+HvUyR7EVcLtD3qAUuIr5h5EjYv99phXta32/+4Eza6fWfXi4X55sU4CLivnPnoG9fqFIFqlUDYPHOxXz353d8VP0jcmTI4XKBvkkBLiLuGz4cDh+GKVPAGGJtLJ1/6Mzdd9xNh4c6uF2dz1KAi4i7zpyBfv2genWnBQ58ufFLNhzewKR6kzRp5wb0JaaIuGvIEPj7b2f4IM4OOD2W9KBcnnI0KtnI5eJ8m1rgInJ7XbgAp0/DqVNw/Dh8/DHUquWMPgE+WfUJB04fYFL9SQQYtTFvRAEuIgl35gwcOABHjzoBfOrUlTA+ffrff776WFRU3PcMCID33wfgyJkj9Pu5H88We5Yqd1dx4QP6FwW4iDiheviwE84HD165v/rnAwecAL6RjBkhUybInPnKfY4ccY9d+3OxYnDffYAzaScyOpL+1fvfhg/t/xTgIqnB2bOwfLkzzvraUD540GlRXytNGsiTx7mVLAk1a0LevM7jXLmcEL46iDNmvLx2ya3QpJ3EU4CLpGTnzsHnn0P//nFDOmdOJ4zz5oXy5Z1QvhTOl+6zZ09SICdW10VdCUkbokk7iZCgADfGZAFGASUBC7QEtgFTgALAbqChtfZ4slQpIokTGQkjRsCHHzpdIzVqQOfOzjKtuXJd3l/SVyzZtYS52+fSr1o/QkNC3S7HbyT0n9dBwAJrbTGgFLAV6AYsttYWBhZ7HouImy5ccCbF3HsvdOzo9C8vXw4LFzpdIPnz+1x4x9pYOi/szF133MWrD73qdjl+5aYBbozJDFQBRgNYa6OstSeAZ4DxntPGA88mV5EichMXLzpriRQpAq+8AgUKwJIlsHQpPOLbS7BO/G0i6w+vp+9jfUmfJr3b5fiVhLTACwERwFhjzHpjzChjTAhwp7X2EIDnPmd8LzbGtDXGhBtjwiMiIrxWuIjg7F4zdiwULQpt20Lu3E5re8UKePRRt6u7oYsxF9l3ch89lvQgLE8YTe5v4nZJfichfeBBwINAB2vtamPMIBLRXWKtHQGMAGdLtVuqUkTiiomBSZOc8dN//eVsgDB0qLMMq0s71sTExnDs/DGOnj1KxLkIjp49evkWcTaCo+eu+vnsUY5HXvnKbGLdiZq0cwsSEuD7gf3W2tWex9NxAvyIMSa3tfaQMSY3EM84JBHxqpgYmDYN3n0Xtm2DUqVg9myoXTtZg9tay8HTB9kSsYUtEVvYcXxHnIA+evYo/5z/h1gbe91rDYbsGbKTMyQnOUNyUipXKUIzhF5+/MCdD1Axf8Vkqz0lu2mAW2sPG2P2GWOKWmu3AdWALZ5bc6Cf5352slYqkprFxsLMmU5wb97sjMueMQOefdarQ/1ibSx7Tuy5HNRb/956+f7UhVOXz8uUNhO5MuYiZ0hOimQvQqX8lS4Hcs6QnISGXAno7OmzExgQ6LUa5YqEjgPvAHxljEkL7ARa4PSfTzXGtAL2Ag2Sp0SRVMxap4X9zjvw22/OMMApU+C555IU3NGx0ew4tiNOSG+J2MIff//B+ejzl8+7M+ROSoSWoNkDzSieozglQktQIrQEOUNyanNhH5CgALfWbgCu29IepzUuIslh3Tpo3RrWr4fCheGrr6BRIwhMXGvWWsvaQ2uZu33u5aDe/s92LsZevHxO/sz5KRFagqoFql4O6uKhxcmWPpu3P5V4kWZiiviiVaucLyQzZ4Zx46BpUwhK3F/Xk5EnmfT7JEasG8GGwxswGAplLUSJ0BI8XeTpy0FdLEcxMqXLlDyfQ5KVAlzE1yxbBk895UxnX7zYmXyTQNZaftn/CyPXjWTK5imcu3iOUneWYmitoTx///NkCc6SjIXL7aYAF/ElixZBnTrORJzFi51x3Qlw7PwxJv42kRFrR7A5YjMZ02ak6f1NafNgG8LyhKm/OoVSgIv4innzoF49ZzblokXOglM3YK1l+Z7ljFw3kulbpnMh5gLl8pRjxNMjaFyysbpFUgEFuIgvmD0bGjSA++93ZlJmz/6vp0acjWD8xvGMWjeKbf9s4450d9D6wda0ebANpXKVuo1Fi9sU4CJumzYNnn8eypaFBQsgy/X91LE2liW7ljBy3UhmbZ3FxdiLVMpfie6Vu9PgvgZkSJPBhcLFbQpwETdNnAjNm0PFivDdd86ok6scOn2IsRvGMnr9aHYe30m29NloV64dbcq2oURoCZeKFl+hABdxy5gxzjjvRx+FOXMgJOTyU9Za+v3Uj7eXvk2MjaFqgar0ebQPdYvXJTgo2MWixZcowEXcMHy4s+zr44/DrFmQ/soyqtZaui3qxkcrP6LhfQ3p82gfCmcv7GKx4qsU4CK326efQqdOzgJUU6dC8JUWdayNpf289gwPH87LYS/zWa3PtEqf/Cv9ZojcTv37O+Fdvz5Mnx4nvKNjo2kxuwXDw4fTpWIXhtYaqvCWG9Jvh8jtYC289x506wZNmsDXX0PatJefjoqJovH0xkzYOIH3q75P/+r9NflGbkpdKCLJzVro0cPZYPjFF2HUqDgLUp2/eJ76U+sz/6/5DKw5kE4Pd3KvVvErCnCR5GStsxv8wIHOlmfDh8dZBvb0hdPUnlyb5XuWM+LpEbQp28bFYsXfKMBFkktsLHToAMOGOfeDBsXZNefY+WM8+dWTrD24lon1JvL8/c+7WKz4IwW4SHKIjYWXXnK6Szp3ho8+ihPeR84coebEmvzx9x/MaDiDZ4o942Kx4q8U4CLgdHV460vD6Gho2RK+/NLp++7dO8577zu5j+pfVmffyX3MbTKXGvfU8M51JdVRgEvqduECtGvnbJqQPj1kzHj9LSQk/uP/9vx77znbnvXuDT17xrncjmM7qDahGscjj7Ow2UIq31XZnc8tKYICXFKviAioWxd+/tlpMd9xB5w5A2fPOvdnzsCxY7B375XHZ85AVNTN3/ujj6BLlziHtkRsofqE6lyIucCSF5ZQNk/ZZPpgkloowCV12rTJmQl5+LDTWm7YMOGvjYq6EvJXh/2lW65cULVqnJesO7SOxyc+TlBAEMteXEbJnCW9+3kkVVKAS+ozbx40bux0fSxbBuXLJ+71adM6t6xZE3T6yn0refKrJ8kSnIVFzRZpXRPxGs3ElNTDWmcdktq14d57Yc2axId3Ii3euZgaX9bgzpA7WdFihcJbvEoBLqnDxYvw3/8665A88wysWAH58iXrJb/d9i1PTXqKQlkLsbzFcu66465kvZ6kPgpwSfmOHYMnnoARI6B7d2cRqavW3k4OUzZNod7Uetx/5/382PxHcmXMlazXk9RJfeCSsm3fDk8/DXv2wPjx8MIL152y8/hO5v05jzQBaQgOCr7pLV1Quis/B6a7btGp0etG0+bbNlS+qzJzn59L5nSZr7umiDcowCXlWrwYnnsO0qSBJUugUqXrTpm6eSqt57TmdNTpW75MusB0ccJ978m91LynJrMazdJelZKsFOCSMn3xhTNBp1gxmDsXChSI83RkdCRvfP8Gw8KHUSFfBcY+M5bM6TITGR1JZHQkF6IvXP75324XYuI/J2+mvPT6Ty/SBaVz57NLqqEAl5QlOtpZe2TQIKhVCyZPvm6j4B3HdtBwekPWHVrHGw+/wYfVPiRNYBqXCha5dQpwSTlOnnQ2S5g/3xlt8vHHcdbdBpi+ZTqt5rQiwAQwu/Fs6hSt41KxIkmnAJeUYedOZ3z39u3OaJM2cdfVvhB9gc4LO/PZms8on7c8U56bQoEsBdypVcRLFODi/1asgHr1ICYGFi6ERx+N8/TO4ztpNL0R4QfD6VShE/2q9yNtYNp/eTMR/6EAF/82frzT2i5Y0PmysnDcmY4zt86k5eyWAMxqNItniz3rRpUiyUITecQ/xcY6GwS/+CJUqQK//BInvKNioug4vyP1p9anSPYirH9pvcJbUhy1wMX/bNrkDBFcvhxeftkZcZLmyiiS3Sd203BaQ9YcXEPHhzryUY2P1GUiKZICXPzHqVPOZgmDBjlrd48eDS1axNntZvYfs3lx9otYa5nRcAb1itdzsWCR5KUuFPF91jrjuYsVg08+gVatnNEmLVteDu+omChe//51np3yLPdkvYd1L61TeEuKpxa4+LYtW6B9e1i6FMqWhW++uW4J2D0n9tBoeiNWH1hN+3LtGVBzgGZBSqqQ4Ba4MSbQGLPeGDPX87igMWa1MeZPY8wUY4w6GcV7Tp+GN9+EUqVgwwYYPhxWr74uvL/d9i1lvijD1r+3Mq3BNIbUGqLwllQjMV0oHYGtVz3uD3xirS0MHAdaebMwSaWshalToXhxZyZl8+awbZuzlvdVsyovxlyky8Iu1Pm6DgWyFGBt27U8V+I5FwsXuf0SFODGmHzAU8Aoz2MDPAZM95wyHtAYLUmaP/6AGjWgUSPImRNWroRRoyA09PIp1lq+3fYtpT4vxYBVA3gl7BVWtlrJvdnudbFwEXcktAX+KfAmEOt5nB04Ya2N9jzeD+SN74XGmLbGmHBjTHhERESSipUU6swZZ0z3Aw/A2rUwdKiz3dnDD8c5LfxgOI+Of5Q6X9chxsYwp/Echj41lOCgYJcKF3HXTb/ENMY8DRy11q41xlS9dDieU218r7fWjgBGAISFhcV7jqRS1sKMGc7CU/v3O5Ny+vd3Wt9X2X1iN28tfovJmyYTmiGUobWG0ubBNlpBUFK9hIxCqQTUMcbUAoKBzDgt8izGmCBPKzwfcDD5ypQUZ/t26NDBWbukVCn4+uvrNlw4fv44fVf0ZfCvgwkwAbxV+S26Vu6qHW5EPG4a4Nba7kB3AE8LvLO1tqkxZhrwHPA10ByYnYx1Skpx9iz07et8QZk+PQwe7MymDLryq3gh+gLD1gyj9/LenIg8QfPSzen9aG/yZU7eTYhF/E1SxoF3Bb42xvQB1gOjvVOSpEiRkU4r+513YO9eZ2/K/v0h15XNfq21TN08le6Lu7PrxC5q3lOTj6p/RKlcpVwsXMR3JSrArbU/Aj96ft4JlL/R+SIcOACff+5scRYR4XSXTJwIjzwS57QVe1bQ+YfO/HrgVx648wG+/7/vqXlPTZeKFvEPmokp3metszrg4MEwfbqzTnft2vDqq/DYY3HWLtn29za6LurK7G2zyZspL2OfGUuzB5oRGBB4gwuICCjAxZsuXIBp05zFpsLDnb0oO3RwpsIXKhTn1KNnj/Luj+8yYu0IMqTJwAePfcBrFV7TLu4iiaAAl6Q7fNjpJvn8czhyBIoWdcZyv/ACZMwY59RzF88xcNVA+v/cn8joSP4b9l96/acXOUNy/subi8i/UYDLrVuzxmltT50KFy/CU0853STVq0NA3Dli5y6eY9yGcXyw4gMOnj5I3WJ16Ve9H0WyF3GpeBH/pwCXxImKcibfDB7s9HNnyuQMA2zf/rrtzMDpKvns188YtmYY/5z/h4r5KzLluSlUvquyC8WLpCwKcEmYI0ec3d6HD4dDh5ywHjzYWWwq8/UTa7b9vY2BqwYyfuN4omKiqFO0Dp0rdqZS/koYE99EXhFJLAW43Nhff0Hv3s4Y7qgoePxxZ4GpJ564rpvEWstPe39iwKoBzNk2h3SB6Xix9It0qtCJojmKuvQBRFIuBbj8uxUr4JlnnOBu08bpJilW7LrTYmJjmPXHLAasHMDqA6vJnj47var0ol35dvpyUiQZKcAlflOmOKNIChaE+fOd+2ucjTrL2A1j+eSXT9h5fCf3ZruXYbWG0bx0cw0HFLkNFOASl7UwYICzG07lyjB7NmTLFueUI2eOOF9Mhg/j2PljPJzvYQbUGECdonU0AUfkNlKAyxXR0dCxIwwb5myqMG4cBF9Za3trxFYGrhrIl799SVRMFHWL1+WNh9+gYv6K7tUskoopwMVx9iw0aQLffuu0vj/8EAICsNayfM9yBqwawNztcwkOCqZlmZZ0qtCJwtmvHzYoIrePAlycIYJPPw3r1jkzKF95BXDW427+TXO+3f4toRlCea/qe7wc9jKhIaE3eUMRuR0U4KndH3/Ak0/C0aPwzTfOolPA+kPrqT+1PvtP7efjGh/Trlw70qdJ73KxInI1BXhqdmmYYJo08OOPUK4cAGPWj+GV714hNCSU5S2WUyFfBXfrFJF4JXRTY0lppkxx1izJmdOZEl+uHOcvnqf1nNa0mtOKR+5+hHVt1ym8RXyYAjy1sdbZzqxxYyhfHlauhIIF2Xl8J5XGVGL0+tH0fKQnC5ouUF+3iI9TF0pqcvUwwYYNYfx4CA5m7va5NJvVDIC5TebyVJGnXC5URBJCLfDU4uxZqFfPCe8uXWDyZGLSpqHnkp7UnlybglkKsq7tOoW3iB9RCzw1uHqY4GefQbt2RJyNoMmMJizetZhWZVrxWa3PCA4Kvvl7iYjPUICndJeGCR45ArNmQZ06/LL/FxpMa0DE2QhG1R5FqwdbuV2liNwCBXhKdvUwwWXLsGFhDP31M17//nXyZc7HqlarKJO7jNtVisgtUh94SnX1MMFVqzhTqjhNZzalw/wOPH7v46xtu1bhLeLnFOApzeHD8PrrV4YJ/vwzf2SO4qFRDzFl8xQ+eOwDZjeeTdb0Wd2uVESSSF0oKcX+/c747hEjnA0YWreGIUOYvnMuLWa3IDgomO//73uqF6rudqUi4iUKcH+3axf06wdjxzqTdF54Abp142KhAnRb1I2BvwykQr4KTH1uKvnvyO92tSLiRQpwf7Vtm7Pk68SJEBjotLi7doW772bfyX08P+Exftr7Ex3Kd2BAzQGkDUzrdsUi4mUKcH/z++/Qt6/zJWVwMLz6KnTuDHnycODUAfrN68DIdSMJDAhkUr1JNLm/idsVi0gyUYD7i7VroU8fZ8nXjBmd1nanTpAzJ/tP7affvPaMXDeSWBtL81LN6fFIDwpmvX4fSxFJORTgvm7lSie458+HLFngnXecVne2bOw7uUBCGM4AAAzwSURBVI8Pv3uF0etHE2tjaVG6Bd0rd1dwi6QSCnBfZK2zPnfv3rB0KeTI4XSbtGsHmTOz9+RePpzbg9HrR2OxtCzdku6PdKdAlgJuVy4it5EC3JdYCwsWOC3ulSshd24YOBDatoWQEPac2MOHc99kzPoxALQs05Lulbtzd5a7XS5cRNygAPcVCxZAz55OX3f+/M7elC1bQnAwu0/s5sNvX2fshrEAtH6wNd0qd+OuO+5yuWgRcZMC3G07djhfRn77LRQqBKNGQbNmkDYtu47vou/CvozbOI4AE0CbB9vQrXI3jecWEUAB7p6zZ51x3AMGOItNffSRs9lC2rTsPL6Tvgv6Mn7jeAJMAC+VfYlulbuRL3M+t6sWER+iAL/drIVp05yx2/v2QdOmTnjnycPO4zv5YP4HjN84nqCAIP5b9r90rdxVwS0i8VKA306bNjlDAJcuhdKlYdIkqFyZvSf38v6c1ozbMI6ggCDalWtH18pdyZMpj9sVi4gPu2mAG2PyAxOAXEAsMMJaO8gYkw2YAhQAdgMNrbXHk69UP3biBLz7rrMbzh13ONuatW3LkfN/8+GC1xgePhxAwS0iiZKQFng08Ia1dp0xJhOw1hjzA/AisNha288Y0w3oBnRNvlL9UGwsjBsH3brB33/DSy9Bnz6cCAlkwLJ3+PSXT4mMjqRF6Ra8/Z+3NapERBLlpgFurT0EHPL8fNoYsxXICzwDVPWcNh74EQX4Fb/+Ch06OPeVKsH333P2viIM+XUI/X/uz4nIEzS6rxHvP/o+RbIXcbtaEfFDidrQwRhTACgDrAbu9IT7pZDP+S+vaWuMCTfGhEdERCStWn9w9Ci0agUPPQR798KXX3Jh6SI+u/gz9wy+h+6Lu1MpfyXWv7Ser5/7WuEtIrcswQFujMkIzABes9aeSujrrLUjrLVh1tqw0NDQW6nRP0RHw6BBUKQITJgAnTsTs3UL40pGU3RoMTrM70DRHEX5qcVPzH1+LqVzlXa7YhHxcwkahWKMSYMT3l9Za2d6Dh8xxuS21h4yxuQGjiZXkT5v6VKnu2TzZqhZE/vpp8y0W3j7q0ps/XsrZXOXZUTtEdQoVANjjNvVikgKkZBRKAYYDWy11g686qk5QHOgn+d+drJU6Mt27XK+oJw6FQoUwM6cycKS6emxtBlrD62leI7izGg4g7rF6iq4RcTrEtICrwQ0A343xmzwHHsLJ7inGmNaAXuBBslToo+5tFLg4MEwZw6kTQvvvcfPTSrz1s/vsXzScgpkKcD4Z8fT9P6mBAYEul2xiKRQCRmF8hPwb83Hat4tx4edO+dsXzZkiDMhJ3t26NqV9Q2r0HPLEOZNeodcGXMxtNZQWj/YWluYiUiy00zMm9m925l4M2oUHD8OpUsTPXok88LuYOTmCcyd/SFZg7PSr1o/OjzUgQxpMrhdsYikEgrw+FzbTWIM1KvHX63rMSZgI+M29uLQrEPkzpibXlV60enhTmQJzuJ21SKSyijArxZPN8n5bm8ws0Y+Ru2ZxY+rmhBoAnmqyFO0LtOaJws/SVCA/ghFxB1KH3C6SYYOhdGjL3eTbPz8XUblOczELSM5sewE92S9h76P9aV56eZaq0REfELqDfB4uklOPlebyc/ew+jTywg/+C7pItJRv0R9WpdpzX8K/IcAk6iJqyIiySr1Bfg13SQ2ezZ+7t6UUcXPM3X3d5zfdp4H7nyAIU8O4fn7nydb+mxuVywiEq/UE+BRUTBypLPT+5EjHC1/HxM+bciogA1sO/YlmfZm4oVSL9D6wdaUzV1WE29ExOel/ACPiYHJk6FXL9i1i+V1HmBwzaLMPraS6BObqZS/Et0e6U6DEg0ISRvidrUiIgmWcgPcWmej4B49YNMm1lQtQo8OZfjh1HpynMtBx4c60qpMK4qHFne7UhGRW5IyA3zZMujeHVatYnPYXbw9oByzzqwhR3QO/lfzf7wc9jLp06R3u0oRkSRJWQG+bh289RZ8/z07i93Ju/0rMPH8ajJdPMH7Vd/ntQqvkSldJrerFBHxipQR4Nu3w9tvw9SpHMifhT59KjAqNpw0F0/RpWIX3qz0JtkzZHe7ShERr/LvAN+/H957D8aO5e+s6ejf62E+S7OemNi1tH2wLT2q9NCkGxFJsfwzwP/5Bz78ED77jFNpYvnkjTD+d8dmzkavplnJZrzzn3comLWg21WKiCQr/wrw06fhk09gwADOXzjD0P+WpV/uHfxzYTX1763P+4++T4nQEm5XKSJyW/hHgF+4AJ9/Dh98QNSxCMa0KE3vew9yMDKcx/M9Tp/H+hCWJ8ztKkVEbiv/CPAnniBm2Y9MblSCd8qkZef5DVQKrcTkatOocncVt6sTEXGFXwT4N20eoWet3Ww+t4Uyd5RhXt2RPHHvE5ruLiKpml8E+JdBm4lOn46ptaZSv0R9rQooIoKfBPjI2iPJnC6zNk8QEbmKXySilnQVEbme+iJERPyUAlxExE8pwEVE/JQCXETETynARUT8lAJcRMRPKcBFRPyUAlxExE8pwEVE/JQCXETETynARUT8lAJcRMRPKcBFRPyUAlxExE8pwEVE/JQCXETETyUpwI0xTxhjthlj/jLGdPNWUSIicnO3HODGmEBgKPAkUAJoYowp4a3CRETkxpKypVp54C9r7U4AY8zXwDPAFm8UdrXXXoMNG7z9riIit0fp0vDpp95/36R0oeQF9l31eL/nWBzGmLbGmHBjTHhEREQSLiciIldLSgvcxHPMXnfA2hHACICwsLDrnk+I5PiXS0TE3yWlBb4fyH/V43zAwaSVIyIiCZWUAF8DFDbGFDTGpAUaA3O8U5aIiNzMLXehWGujjTHtge+BQGCMtXaz1yoTEZEbSkofONbaecA8L9UiIiKJoJmYIiJ+SgEuIuKnFOAiIn5KAS4i4qeMtbc0t+bWLmZMBLDnFl+eA/jbi+V4i+pKHNWVOKorcVJqXXdba0OvPXhbAzwpjDHh1towt+u4lupKHNWVOKorcVJbXepCERHxUwpwERE/5U8BPsLtAv6F6koc1ZU4qitxUlVdftMHLiIicflTC1xERK6iABcR8VM+H+DGmDHGmKPGmE1u13I1Y0x+Y8xSY8xWY8xmY0xHt2sCMMYEG2N+NcZs9NT1nts1XWKMCTTGrDfGzHW7lqsZY3YbY343xmwwxoS7Xc8lxpgsxpjpxpg/PL9nD/tATUU9f06XbqeMMa+5XReAMaaT53d+kzFmsjEm2O2aAIwxHT01bfb2n5XP94EbY6oAZ4AJ1tqSbtdziTEmN5DbWrvOGJMJWAs8a631+p6giazLACHW2jPGmDTAT0BHa+0vbtYFYIx5HQgDMltrn3a7nkuMMbuBMGutT00AMcaMB1ZYa0d51tzPYK094XZdl3g2Nj8APGStvdUJet6qJS/O73oJa+15Y8xUYJ61dpzLdZUEvsbZQzgKWAC8bK390xvv7/MtcGvtcuCY23Vcy1p7yFq7zvPzaWAr8ewJertZxxnPwzSem+v/Shtj8gFPAaPcrsUfGGMyA1WA0QDW2ihfCm+PasAOt8P7KkFAemNMEJAB39ghrDjwi7X2nLU2GlgG1PXWm/t8gPsDY0wBoAyw2t1KHJ6uig3AUeAHa60v1PUp8CYQ63Yh8bDAQmPMWmNMW7eL8SgERABjPd1Oo4wxIW4XdY3GwGS3iwCw1h4ABgB7gUPASWvtQnerAmATUMUYk90YkwGoRdytKJNEAZ5ExpiMwAzgNWvtKbfrAbDWxlhrS+PsU1re89841xhjngaOWmvXulnHDVSy1j4IPAm083TbuS0IeBAYbq0tA5wFurlb0hWeLp06wDS3awEwxmQFngEKAnmAEGPM/7lbFVhrtwL9gR9wuk82AtHeen8FeBJ4+phnAF9Za2e6Xc+1PP/l/hF4wuVSKgF1PH3NXwOPGWMmulvSFdbag577o8AsnP5Kt+0H9l/1v6fpOIHuK54E1llrj7hdiEd1YJe1NsJaexGYCVR0uSYArLWjrbUPWmur4HQHe6X/GxTgt8zzZeFoYKu1dqDb9VxijAk1xmTx/Jwe5xf7DzdrstZ2t9bms9YWwPlv9xJrreutIwBjTIjnS2g8XRQ1cf7b6ypr7WFgnzGmqOdQNcDVL8iv0QQf6T7x2AtUMMZk8PzdrIbzvZTrjDE5Pfd3AfXw4p9bkvbEvB2MMZOBqkAOY8x+4B1r7Wh3qwKcVmUz4HdPfzPAW559Qt2UGxjvGSEQAEy11vrUsD0fcycwy/k7TxAwyVq7wN2SLusAfOXprtgJtHC5HgA8fbk1gJfcruUSa+1qY8x0YB1OF8V6fGda/QxjTHbgItDOWnvcW2/s88MIRUQkfupCERHxUwpwERE/pQAXEfFTCnARET+lABcR8VMKcBERP6UAFxHxU/8Pn+BlI9ZdWvMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "iter = 100\n",
    "range_x1 = np.linspace(1,9,17)\n",
    "results = []\n",
    "for sigma in range_x1:\n",
    "    calu = np.zeros(3)\n",
    "    for i in range(iter):\n",
    "        gussian_noise = np.random.normal(0,sigma,size=(num,1))\n",
    "        dectected_slope = slope_distance + gussian_noise\n",
    "        sol4_root = root(f5,[x_axis,y_axis,z_axis])#,method='hybr')\n",
    "        sol4_fsolve = fsolve(f5,[x_axis,y_axis,z_axis])\n",
    "        calu += (sol4_fsolve-gt)**2\n",
    "    calu /= iter\n",
    "    results.append(np.sqrt(calu))\n",
    "results = np.array(results).reshape(-1,3)\n",
    "plt.figure(1)\n",
    "plt.title(\"%d\" %(num))\n",
    "plt.plot(range_x1,results[:,0],'r',range_x1,results[:,1],'b',range_x1,results[:,2],'g')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.  1.5 2.  2.5 3.  3.5 4.  4.5 5.  5.5 6.  6.5 7.  7.5 8.  8.5 9. ]\n",
      "[[10.11902189  0.          9.42666111]\n",
      " [14.49332217  0.         13.7541659 ]\n",
      " [18.03956248  0.         16.98406218]\n",
      " [22.92453384  0.         21.31950166]\n",
      " [26.94352029  0.         25.4359226 ]\n",
      " [34.18724255  0.         31.8402409 ]\n",
      " [37.75055436  0.         35.59649382]\n",
      " [42.73360502  0.         40.15570045]\n",
      " [49.67562627  0.         46.38332158]\n",
      " [49.29132894  0.         46.12750264]\n",
      " [57.32803022  0.         53.99289387]\n",
      " [61.53511121  0.         57.26745338]\n",
      " [62.3907509   0.         58.61508609]\n",
      " [62.82262767  0.         59.31728498]\n",
      " [79.01328043  0.         74.22957848]\n",
      " [73.78313403  0.         69.40305974]\n",
      " [86.00698816  0.         80.36091694]]\n"
     ]
    }
   ],
   "source": [
    "print(range_x1)\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-3-b5958880b1ec>, line 2)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-3-b5958880b1ec>\"\u001b[1;36m, line \u001b[1;32m2\u001b[0m\n\u001b[1;33m    while cmt--:\u001b[0m\n\u001b[1;37m               ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "cmt = 10\n",
    "while cmt--:\n",
    "    print(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[31.48469453  0.39942985 33.33081942]\n",
      "[22.62979084  0.34973027 23.87993926]\n",
      "[1.91635543e+01 4.13041537e-06 2.04958750e+01]\n",
      "[1.56530201e+01 5.49719986e-03 1.65614459e+01]\n",
      "[1.49855590e+01 1.45022016e-05 1.59310082e+01]\n",
      "[1.37735743e+01 7.19550627e-04 1.44625340e+01]\n",
      "[1.23773873e+01 1.34166512e-05 1.32827954e+01]\n",
      "[11.61100686  0.02121365 12.19345787]\n",
      "[1.08648467e+01 4.21813874e-05 1.16488048e+01]\n",
      "[9.86497438e+00 9.89703124e-06 1.05320623e+01]\n"
     ]
    }
   ],
   "source": [
    "count =100\n",
    "results=[]\n",
    "for i in range(10):\n",
    "    results.append([])\n",
    "while (count):\n",
    "    myint = np.random.randint(0,10000)\n",
    "    #print(myint)\n",
    "    for k in range (3,13):\n",
    "        num = k\n",
    "        np.random.seed(myint)\n",
    "        #nois_dist = np.random.\n",
    "        #axis = np.random.uniform(low,high,(num,3))\n",
    "        axis = np.zeros((num,3))\n",
    "        y_axis = np.random.uniform(-20000,20000,size=(1,))\n",
    "        z_axis = np.random.uniform(6000,18000,size =(1,))\n",
    "        x_axis = np.random.uniform(-10000,10000,size =(1,))\n",
    "       # print(x_axis)\n",
    "       # print(y_axis)\n",
    "       # print(z_axis)\n",
    "        axis[:,1] = y_axis\n",
    "        axis[:,2] = np.random.uniform(0,4000,size =(num,))\n",
    "        axis[:,0] = np.random.uniform(x_axis+10000,x_axis+12500,size =(num,))#np.random.uniform(4000,6500,size =(num,)) #np.random.uniform(6000,8500,size =(num,))\n",
    "        #np.random.uniform(8000,10500,size =(num,)) #np.random.uniform(-1400,1100,size =(num,))\n",
    "        #axis[:,0] = x_axis + np.abs(z_axis-axis[:,2]) +np.random.uniform(-10,10,(num,))\n",
    "        #np.random.uniform((z_axis-axis[:,2]) /3**0.5,(z_axis-axis[:,2]) * 3**0.5)\n",
    "        #np.random.uniform((x_axis-1,x_axis+1),3)\n",
    "        gt = np.array((x_axis,y_axis,z_axis)).squeeze()\n",
    "        #slope_distance= np.random.randint(10,20,cal).reshape(cal,1)\n",
    "        slope_distance = np.sqrt((x_axis-axis[:,0])**2 + (z_axis-axis[:,2])**2) .reshape(num,1)\n",
    "        iter = 10\n",
    "        #range_x1 = np.linspace(1,9,17)\n",
    "        \n",
    "        #range_x1:\n",
    "        calu = np.zeros(3)\n",
    "        sigma =3.0\n",
    "        for i in range(iter):\n",
    "            gussian_noise = np.random.normal(0,sigma,size=(num,1))\n",
    "            dectected_slope = slope_distance + gussian_noise\n",
    "            sol4_root = root(f5,[x_axis,y_axis,z_axis])#,method='hybr')\n",
    "            sol4_fsolve = fsolve(f5,[x_axis,y_axis,z_axis])\n",
    "            calu += (sol4_fsolve-gt)**2\n",
    "        calu /= iter\n",
    "        results[k-3].append(np.sqrt(calu))\n",
    "    count -= 1\n",
    "for j in range(10):\n",
    "    results[j] = np.array(results[j]).reshape(-1,3)\n",
    "    results[j] = np.mean(results[j],axis = 0)\n",
    "    print(results[j])\n",
    "        #plt.figure(k)\n",
    "        #plt.title(\"%d\" %(num))\n",
    "        #plt.plot(range_x1,results[:,0],'r',range_x1,results[:,1],'b',range_x1,results[:,2],'g')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0. , 1.5, 3. , 4.5])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linspace(0,4.5,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([ 7.93905316,  0.        , 13.36177116]),\n",
       " array([10.56601648,  0.        , 17.66983455]),\n",
       " array([ 7.4386378,  0.       , 12.5274809]),\n",
       " array([ 9.76486018,  0.        , 16.43237519]),\n",
       " array([5.87627917, 0.        , 9.96960922]),\n",
       " array([ 6.49009916,  0.        , 11.66786325]),\n",
       " array([ 6.88698535,  0.        , 12.12525984]),\n",
       " array([ 5.77138035,  0.        , 10.29004604]),\n",
       " array([5.42950936, 0.        , 9.94326416]),\n",
       " [array([5.29185822, 0.        , 9.19089388])]]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(1)\n",
    "    #nois_dist = np.random.\n",
    "    #axis = np.random.uniform(low,high,(num,3))\n",
    "print(np.random.uniform(-20000,20000,size=(1,)))\n",
    "print(np.random.uniform(6000,18000,size =(1,)))\n",
    "print(np.random.uniform(-10000,10000,size =(1,)))\n",
    "np.random.uniform(0,4000,size =(12,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[50.18070507  1.48067788 55.46994469]\n",
      "[50.18070507  1.48067788 55.46994469]\n",
      "[50.18070507  1.48067788 55.46994469]\n",
      "[50.18070507  1.48067788 55.46994469]\n"
     ]
    }
   ],
   "source": [
    "#####画表版\n",
    "count =100\n",
    "results=[]\n",
    "for i in range(4):\n",
    "    results.append([])\n",
    "while (count):\n",
    "    myint = np.random.randint(0,10000)\n",
    "    #print(myint)\n",
    "    #for k in range (1,5):\n",
    "        num = 12\n",
    "        np.random.seed(myint)\n",
    "        #nois_dist = np.random.\n",
    "        #axis = np.random.uniform(low,high,(num,3))\n",
    "        axis = np.zeros((num,3))\n",
    "        y_axis = np.random.uniform(-20000,20000,size=(1,)) #+1500*k\n",
    "        z_axis = np.random.uniform(8000,18000,size =(1,))\n",
    "        x_axis = np.random.uniform(-10000,10000,size =(1,))\n",
    "        gt = np.repeat(np.array((x_axis,y_axis,z_axis)).reshape(1,3),4,axis = 0)\n",
    "        gt[:,1] += np.linspace(0,4.5,4)\n",
    "        \n",
    "       # print(x_axis)\n",
    "       # print(y_axis)\n",
    "       # print(z_axis)\n",
    "        axis[:,1] = y_axis\n",
    "        axis[:,2] = np.random.uniform(0,2000,size =(num,))\n",
    "        axis[:,0] = np.random.uniform(x_axis+10000,x_axis+12500,size =(num,))#np.random.uniform(4000,6500,size =(num,)) #np.random.uniform(6000,8500,size =(num,))\n",
    "        #np.random.uniform(8000,10500,size =(num,)) #np.random.uniform(-1400,1100,size =(num,))\n",
    "        #axis[:,0] = x_axis + np.abs(z_axis-axis[:,2]) +np.random.uniform(-10,10,(num,))\n",
    "        #np.random.uniform((z_axis-axis[:,2]) /3**0.5,(z_axis-axis[:,2]) * 3**0.5)\n",
    "        #np.random.uniform((x_axis-1,x_axis+1),3)\n",
    "        \n",
    "        #slope_distance= np.random.randint(10,20,cal).reshape(cal,1)\n",
    "        slope_distance = np.sqrt((x_axis-axis[:,0])**2 + (z_axis-axis[:,2])**2) .reshape(num,1)\n",
    "        iter = 10\n",
    "        #range_x1 = np.linspace(1,9,17)\n",
    "        \n",
    "        #range_x1:\n",
    "        calu = np.zeros(3)\n",
    "        sigma =3.0\n",
    "        for i in range(iter):\n",
    "            gussian_noise = np.random.normal(0,sigma,size=(num,1))\n",
    "            dectected_slope = slope_distance + gussian_noise\n",
    "            sol4_root = root(f5,[x_axis,y_axis,z_axis])#,method='hybr')\n",
    "            sol4_fsolve = fsolve(f5,[x_axis,y_axis,z_axis])\n",
    "            calu += (sol4_fsolve-gt)**2\n",
    "        calu /= iter\n",
    "        results[k-1].append(np.sqrt(calu))\n",
    "    count -= 1\n",
    "for j in range(4):\n",
    "    results[j] = np.array(results[j]).reshape(-1,3)\n",
    "    results[j] = np.mean(results[j],axis = 0)\n",
    "    print(results[j])\n",
    "        #plt.figure(k)\n",
    "        #plt.title(\"%d\" %(num))\n",
    "        #plt.plot(range_x1,results[:,0],'r',range_x1,results[:,1],'b',range_x1,results[:,2],'g')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x20b79b728c8>,\n",
       " <matplotlib.lines.Line2D at 0x20b79bad8c8>,\n",
       " <matplotlib.lines.Line2D at 0x20b79badd48>]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hUZfrG8e+ThARCgCBFkCKirivwU9BgFxCsiGJBVFSUlcUuIq6KroroWrCBFViqDVSwoCICojSlBASkKotKVSLSQkuZ5/dHxl2MQRKY5GQm9+e6cjEz582cW4SbN++cYu6OiIhEv7igA4iISGSo0EVEYoQKXUQkRqjQRURihApdRCRGqNBFRGKECl1EJEao0EXyMbMjzWyXmb0edBaRolChi/zRS8CcoEOIFJUKXWQPZnYFsBn4LOgsIkWlQhcJM7PKQB+gZ9BZRPaHCl3kfx4Bhrj76qCDiOyPhKADiJQGZtYUOBNoFnQWkf2lQhfJ0wpoAKwyM4AUIN7MGrn7cQHmEik00+VzRcDMkoHKe7x0F3kFf5O7ZwQSSqSINEMXAdx9B7Djt+dmlgnsUplLNNEMXUQkRugoFxGRGKFCFxGJESp0EZEYoUIXEYkRgR3lUr16dW/QoEFQuxcRiUpz5879xd1rFLQtsEJv0KAB6enpQe1eRCQqmdmPe9umJRcRkRihQhcRiREqdBGRGKFCFxGJESp0EZEYoUIXEYkRKnQRkRihQhcRKUF9Hm7NzPGDi+W9VegiIiXky3EDeYjPGT/z9WJ5fxW6iEgJCOXm0GNiT2pvj+Mf3d8uln3ojkUiIiXgrcE9mJ26nWHVrqdi1ZrFsg/N0EVEitnOrb9y74pXaLq5Ap1vGlBs+9EMXUSkmPV/vhOrUnIZ9n+PEBdffLWrQhcRKUYbfljMYzs+5cKdtWh9Sc9i3ZeWXEREitFDr3RkZwL07TS02PelQhcRKSaLv/yAQeWXcFPWMRzV/Lxi358KXUSkmNz1TlcqZxkP3Tq6RPanNXQRkWLw6Vv/YnzqLzyd2I5qdY8skX1qhi4iEmE5WbvoOesRGm5L4Nbb3yix/WqGLiISYUNf7sriKrsZXe8ukipWLrH9aoYuIhJB2zau44H1b3Lapspcct2TJbpvzdBFRCLoif6XsSHZ+bBNPyyuZOfMKnQRkQhZteQrng19yVXbG3DC2V1KfP9achERiZD7hnQC4LHr3wxk//ssdDMrb2azzWyBmS02s4cLGJNkZm+Z2Qozm2VmDYojrIhIaTV7wjDeqPwDd8adQv1GJweSoTAz9N1Aa3c/FmgKnGtmJ+Ubcz2wyd2PAJ4DSvaTABGRAHkoxJ3julNzh3Fv93cCy7HPQvc8meGn5cJfnm9Ye2BE+PFooI2ZWcRSioiUYu8Ov4cZVbfxSO1OVKp2SGA5CrWGbmbxZjYf2ABMdPdZ+YbUAVYDuHsOsAWoFsmgIiKl0e7tW7l7ST+abEnibzcXz71CC6tQhe7uue7eFKgLnGBmTfINKWg2nn8Wj5l1M7N0M0vPyMgoeloRkVLmxeevYmWlHJ456SESEssHmqVIR7m4+2bgC+DcfJvWAPUAzCwBqAL8WsD3D3L3NHdPq1Gjxn4FFhEpLX5ZvZxHtn7EuZurc3bHXkHHKdRRLjXMLDX8uAJwJrAs37CxwLXhxx2Aye7+hxm6iEgs6fNSR7YlwtOXBbvU8pvCnFhUGxhhZvHk/QPwtrt/ZGZ9gHR3HwsMAV4zsxXkzcyvKLbEIiKlwLLZ43g5aSHddjWi8Sntg44DFKLQ3X0h0KyA1x/c4/Eu4LLIRhMRKb3uHvk3kivAwze9HXSU/9Kp/yIiRTT53Wf4MPVnHo8/h5oNGgcd57906r+ISBHkZmfRc9oDHLotnjvuGBV0nN/RDF1EpAheHXAj81N3MrL2bZRPSQ06zu9ohi4iUkjbN23g/tUjOHFzRS7v2i/oOH+gGbqISCE91b8j6yuGGN38mRK/1nlhqNBFRAph7bfp9M2ZQscd9Til7Q1BxylQ6fsnRkSkFPrnoCvJNXjiuteDjrJXKnQRkX2YN/kNRqSsoHuoOYcd0yLoOHulQhcR+RM/LJrOLR/eSLVdxn23B3et88LQGrqISAHmTX6Dpz6+j3dSVmGVYGitG0g9+NCgY/0pFbqISJiHQkx453Ge+uoZPqu6iUpJ0CMnje6dX6buUc2DjrdPKnQRKfOyd+3grWE9efrb4SxI3cUhiXE8mXAeN9w0gCo16wcdr9BU6CJSZm3buI7B/76Z5zZ+xOqUXBpZEkMP6kKnu/qRVLFy0PGKTIUuImXOTysX8vywG3klZyabyzsts6vwyhF3cN4V/yQuPnprMXqTi4gU0bLZ43jmnTt5tfxysuPh0u11+McZj3DC2V2CjhYRKnQRiXkzPn6FvpP6MDb1J8onwfVZjbnziv4c0axN0NEiSoUuIjEplJvDB6/dz1PzX+arqplUSzIe9Bbc2m0ANeofHXS8YqFCF5GYkz5xBNd+cgNLquzmsIQEXkjuQJfbX6Ji1ZpBRytWKnQRiRm52Vn0feICHsyeQK34eEYdcjuXXvskCYnlg45WIlToIhITVi3+kmv+fR5Tq26lY2Y9BvxjClVrHxZ0rBKla7mISNQbNfA2jnn9VOZV3MqI6n9n1NM/lLkyh0IUupnVM7PPzWypmS02s+4FjGllZlvMbH7468HiiSsi8j9bN6ymc8+GXPnTizTamcKCjp/T+ZZBpfLmEyWhMEsuOUBPd59nZpWAuWY20d2X5Bs3zd3bRT6iiMgfzfjoZa7+4nZWp+TSm1bc3/eTMrNWvjf7LHR3Xw+sDz/eZmZLgTpA/kIXESl22bt28Mjj5/AvptOABKadNJCTz+sWdKxSoUg/l5hZA6AZMKuAzSeb2QIz+8TMGu/l+7uZWbqZpWdkZBQ5rIiULts2ruORPm34d7/OZKxaWuz7W/H1Z5x+b00eiZtO5+1HML/XjyrzPZi7F26gWQowBfiXu7+bb1tlIOTumWbWFujv7kf+2fulpaV5enr6fsYWkaCtXPAF7Uecy6IquwGIC0HLLal0qHcOF190L7UPbxqxfXkoxPAXu3Lbz8MoFzIGHt6djl2fi9j7RxMzm+vuaQVtK9QM3czKAWOAN/KXOYC7b3X3zPDjcUA5M6t+AJlFpBSb/O4zNB/ZmrVJWUxo9DgLWr/D/ZzO+rgd3JL5FnVea0aLO6rQv+8lrF42+4D29evaFVx2V33+tmkYzbensrDzV2W2zPdlnzN0MzNgBPCru9+xlzG1gJ/d3c3sBGA0cKj/yZtrhi4SfTwU4sVnOtIjcwx/3ZbEB53HcXjT1r8bs2TmWEZ/2o/RW77imyq7ADhpcwodarbi0gvuoUGT0wq9v8ljnqbzzHvYUCHEo0nn0fPu94kvlxjR/6Zo82cz9MIU+mnANOAbIBR++T6gPoC7DzCzW4GbyDsiZidwp7t/+Wfvq0IXiS67t2/l5t5pDE35jvaba/HaP+dSqdohf/o936Z/ypjxzzJ64zTmpe4E4PjNyXSo3oIObe/a68Wxdm/fygOPncnT5ebwl22JvHn+UI5rfVXE/5ui0QEVenFRoYtEj59WLuSS50/lq6qZPBA6nd4PTi7ydcNXLviCMR8/xeiMKcxO3Q7AsZvL06HqqVx6bg+OPuF8AJbO+ohOozoyP3UnN+5oxDP3TSG5ilZwf6NCF5H9lj7pVS76tAubEkMMb9iDy65/9oDfc9WSr3j3w76MXv8ZM6puA6DRliROT2jIiMSlpOQYQxvfzwXXPHLA+4o1KnQR2S9vDLiZrmte4eBd8bx/wRs0bXl5xPex9tt03vvgScasncTUKps5e2t1ht32GbUaHhPxfcUCFbqIFEludha9ep/GU4lzaLmpCu/c+VWJXEN89/atUXkvz5L0Z4Wuqy2KyO9s/vlHOj3RnE9SM7h5RxP6PTGLcuWTS2TfKvMDUzavYCMiBVo2exwn9j2SiZUyGFjlKl568psSK3M5cJqhiwgA4958mCsX9SYpwZjc/AVOv+DWoCNJEanQRco4D4Xo+/j59MoeT9OdFXj/759Rv9HJQceS/aBCFynDdmz5ha590hhZ+Ucu31qPob3n6ZjvKKY1dJEyavXSWZzeuz6jKv3IY/FnM/LpH1TmUU4zdJEyaMbHr3Dp1FvZUSHE2CMeoN3VfYKOJBGgQhcpQ9Ysn8PDg69hWPJyDsspx+eXfMDRJ+pGY7FChS5SBvyyejlPDLiKF+Pm4uXhlqxj6X3Xe2XyRsqxTIUuEsO2bVzHcy9ezdO7P2d7Obhm++H0vm54kS5hK9FDhS4Sg3Zv38qAl67jXxvfJyPZuXhnbR69cACNTrow6GhSjFToIjEkNzuL1wbezEM/DmdVSi5n7E7l8dZPc+I51wcdTUqACl0kBngoxPuv9uL+hf1YWiWL43OSGdzoAc689G4sTkcnlxUqdJEoN/ndZ+g17SFmp27nryQyut5dXPLAkyryMkiFLhKl5kwcwX0f38mkqr9SLyGeIVWvo3OvV0hILB90NAmICl0kyiyd9REPvHUjY6qspXp547nyF3Fjz2GUT0kNOpoETIUuEiV+XDyDh4d3YUTydyRXgN60okfP16hco27Q0aSUUKGLlDLbN21g6dzxLF42lUXrF7B4x48sSviV1Sm5JJaH7tnH0evG10vkDkISXfZZ6GZWD3gVqAWEgEHu3j/fGAP6A22BHcB17j4v8nFFYseuzM0sn/spi5ZOYfHa+SzKXMni+I18n5KDW96YpDg4mgq0CNWjSfxRdLrsIV3aVvaqMDP0HKCnu88zs0rAXDOb6O5L9hhzHnBk+OtE4JXwryJlXvbO7Xz79UQWL/6CRWvmsXjrf1hsGXxXKZtQ+ECUBIOjLIk0r811cUfSuN7xNPm/NjQ8pqU+5JRC22ehu/t6YH348TYzWwrUAfYs9PbAq553x+mZZpZqZrXD3ytSJmTt2MZ38z9jyZIpLFnzNUu2/ofFZPBtym6y4/PGxAFHxJWjsdegox1O4zrH0aTRGRzZrA2JFVICzS/Rr0hr6GbWAGgGzMq3qQ6weo/na8Kv/a7Qzawb0A2gfv36RUsqUkrsytzMt/MmsmTpVJasm8/ibStZwi98VymL3PCM24CGceU4OvcgLvCGNDm4KY2Pbslfjz9HR6NIsSl0oZtZCjAGuMPdt+bfXMC3+B9ecB8EDAJIS0v7w3aR0mTHll9Y/vVEliybxpK1C1icuZIlcRv5T8r/lkrigCOsHI28Opf4YTQ+uCmNjm7BUc3OokLlgwLNL2VPoQrdzMqRV+ZvuPu7BQxZA9Tb43ldYN2BxxMpWT9/v4juL7RlTtz63304mWBwpCVyjNfkSjucRofkFfdfmp1FUsXKwYYWCSvMUS4GDAGWuvuzexk2FrjVzEaR92HoFq2fS7RZ+206bQaeyuoKWbTbWY9r7Qga1W1Go8YtOeLY1lrjllKvMDP0U4FrgG/MbH74tfuA+gDuPgAYR94hiyvIO2yxS+SjihSfVUu+ovWQFmxIyuHTk1/itHY3Bx1JpMgKc5TLdApeI99zjAO3RCqUSEn6fsEUWr/Whk2JuUxsOViXmpWopTNFpUz7bu4EWo9qy46EEJ+d+RrHt7k66Egi+02FLmXW0pkf0vrdi8iNcya3HcWxLToGHUnkgKjQpUxaOG00Z37UkXiMLy5+X7dmk5igQpcyZ97kNzhrwjVUCBmTLx/HX9LOCTqSSEToliZSpsz6dAhtJl5NpZw4pl79mcpcYooKXcqM6R+9zFlTunJQdgJTukyh4bGtgo4kElEqdCkTPn/vWc796hZq7y7H1G5fcWjjU4OOJBJxKnSJeRPeeoy2c3ty6M4kptySTp2/pAUdSaRYqNAlpn30+oNcsOh+jtpeni/u+JpaDY8JOpJIsdFRLhKz3ht2D5d/35djMpOZ8I+FHHTI4UFHEilWKnSJSW8N6s5Va5+n+bYUxvdaTJWauv6+xD4VusScV1/qRpcN/+bULZX5+IGlVKp2SNCRREqECl1iyuD+19Jt06ucsaUqYx9aRsWqNYOOJFJiVOgSM156uiO3bn+Hc7dU590+y3XHIClzVOgS9bJ2ZtL/ucu5O3scF2w+mHce/VZ3EZIySYUuUWdrxhq+mjyC6UvGM23LN8xK2cKucnDpljq8+dgy3VlIyiwVupR6P61cyPTPRzBtxWdM3/Ut8yvvJBSXdxJFM0vmxtxmtDjiLC7o9DAJieWDjisSGBW6lCoeCrHi68+YNuNNpq2aznT/kRWVswGoEA8neSr3WxqnN2rLSa076wgWkT2o0CVQOVm7WDBtNNPSxzD95zlML7een5NDAByUYJy2+2BuSEzj9OMuolnLK0gsXzHgxCKllwpdAvHTyoX8/YWz+CJ5A5mJea81iEvgrNxDOb3KqZx2Ukf+ekJb4uLigw0qEkVU6FLidm/fyqXPn8bXKdvoktOE02q34vSWnal7VPOgo4lEtX0WupkNBdoBG9y9SQHbWwEfAN+HX3rX3ftEMqTEDg+FuO3hE/my6jbeqnMHHbs+F3QkkZhRmBn6cOBF4NU/GTPN3dtFJJHEtAHPXcW/Ky6jV+4pKnORCNvn5XPdfSrwawlkkRg3dewL3L5lFOdvrskjD34RdByRmBOp66GfbGYLzOwTM2u8t0Fm1s3M0s0sPSMjI0K7lmjw4+IZdJjRncO3J/JGrznEJ5QLOpJIzIlEoc8DDnX3Y4EXgPf3NtDdB7l7mrun1ahRIwK7lmiwY8svXDz4LHbHOx9c8b4uZStSTA640N19q7tnhh+PA8qZWfUDTiYxwUMhru9zPPOr7GRkk94c1fy8oCOJxKwDLnQzq2VmFn58Qvg9Nx7o+0ps6Pv4+YyqvIrHEs6hbaeHgo4jEtMKc9jiSKAVUN3M1gAPAeUA3H0A0AG4ycxygJ3AFe7uxZZYosa4Nx+mV/Z4Lt9aj3ueHhd0HJGYt89Cd/cr97H9RfIOaxT5r2/TP6XTN705dlcFhjw0F4vT/chFipv+lknEbdmwivZvXkg5N97vOpGKqfoAXKQk6NR/iahQbg5XP96cFZWymHR8Pw5tfGrQkUTKDBW6RNSDD7fio9QNvJjSkZbtuwcdR6RM0ZKLRMw7Q+7kX/Ez6Jp5FDf3GBl0HJEyRzN0iYgFU9/mupXPcfL2FF7812x9CCoSABW6HLBfVi+n/dhOVLU4xtw2TTdoFgmICl0OSPauHXR85iR+qpTLtBbDqX1406AjiZRZKnQ5ID17n8znVTczovrfaX7WtUHHESnTtNAp+23o8114ocJCeuw+js63DAo6jkiZpxm67JeZ4wdzU8Zwztx2EH37zgg6joigQpf9sO67eVwy+QbqhhJ46+7ZJCSWDzqSiKBClyLalbmZi19uwdbkEBPajuagQw4POpKIhKnQpdA8FOKm3s2Znbqddw+9myanXhx0JBHZgz4UlULZtnEdnf5xGMMrreBBb8HF1z0ZdCQRyUczdNmnRTPeo8O7V/BdShb/ij+be3t9HHQkESmACl3+1Gsv38gN6wZSOSGOSU2f5YyLewQdSUT2QoUuBdqVuZnb+5zMvysuo8X2Koy69QudBSpSyqnQ5Q/+M38yHUa0Y37qTu7JOYlHn/xchyaKRAEVuvzOe8Pvocvyvlh548PDH6Dd1X2CjiQihaRCFyDvIlv39mnBs0lzSduVzNt/G89h/3d60LFEpAj2ediimQ01sw1mtmgv283MnjezFWa20MyOi3xMKU5rv03njHtr8WzSXG7e0YTpj65XmYtEocIchz4cOPdPtp8HHBn+6ga8cuCxpKRMfOcJmg49gfkVt/Fm7Vt56clvdD1zkSi1zyUXd59qZg3+ZEh74FV3d2CmmaWaWW13Xx+hjFIMcrOzePSxc3jYv+DorCRGXz6ao09sF3QsETkAkVhDrwOs3uP5mvBrfyh0M+tG3iye+vXrR2DXsj8yVi3l6mdPZ0LVjVydeRgDHphJxao1g44lIgcoEqf+WwGveUED3X2Qu6e5e1qNGjUisGspqi/HDaTZC02YUmkjA6tczat9V6jMRWJEJGboa4B6ezyvC6yLwPtKBHkoRL++l3D3jg+o5wl8ecarHNf6qqBjiUgERWKGPhboHD7a5SRgi9bPS5ctG1bR4a563Ln7A87PrMW8f6xQmYvEoH3O0M1sJNAKqG5ma4CHgHIA7j4AGAe0BVYAO4AuxRVWCi8naxfpk19n0uxRDNs6hR8r5fBUYjt6Pv0BFqeLbIrEosIc5XLlPrY7cEvEEsl+8VCIZXPGMWnqcD5bN4PPK/zE1qS8bcfnJjOieX9Oa3dzsCFFpFjpTNEotu67eXw2cRCTVk5kUtwPrKsYAuCw+AQuz/krZzY8lzPO7kaN+kcHnFRESoIKPYpszVjDlPEDmbTkQybtXsaSKrsBqJZgtMmqw5lVWtHmjOtpeGyrYIOKSCBU6KVY1s5MZk4cxmfzRjNp89fMqrKN3DioEA+nh6pxXbk2nHnK1Rx7+mXExet/pUhZpxYoZXKyd/PB6/9k6MIRTEnOYHsixDk0pyL3+imceUwHTj7rbySlVAk6qoiUMir0UmLT+u8ZPOw2Xtw0nlUpudRPiOe63P/jzAbn0+rcG0k9+NCgI4pIKadCD9iy2eN4fszdjEhYzI5EaJldhX71b+DCq/oQn5gUdDwRiSIq9ACEcnP49O3H6D/7eT5N3UhSInTaeQS3n/0oTVteHnQ8EYlSKvQSlPnrT4wY1p0X1r3P8spZ1CoXxyPWmm7XP0/NBo2DjiciUU6FXgK+XziVF9+6kyGhuWwpD81DFXn94Ou57Nq+JFZICTqeiMQIFXox8VCIqR++QL8vnmRs5fVYAnTYVo/uLf7JSed01en3IhJxKvQI25W5mZHD7qT/ypEsSN3FQUnGPaGTufnq/tQ9qnnQ8UQkhqnQI2T10lkMGnUXA3fNICPZaWxJDEq9hqt6PEdy5WpBxxORMkCFfgBys7P45K1HGTh3IOMqb8AN2mUdTPcT7qL1xXdqWUVESpQKfT+sWT6HIaPuZnDmNNak5FKrXBz3+il0vfRxDjumRdDxRKSMUqEXUm52FuPffoyB6QP4uPLPhOLg7Jxq9K//Ny64sjflyicHHVFEyjgV+j6sWT6HoaPuYXDmVFan5HJwYhz3+Cl0vfgxGh7TMuh4IiL/pUIvQG52Fp++8zgD5wzgo8o/EYqDs3Kq8Vy9LlzY6WHNxkWkVFKh72Htt+l5s/FtU1iVkkvNxDjuDp3M3y95TNcYF5FSr8wXem52FhNGP8HA2a/wUeWfyI2Ds7IP4pm613Fhp4d1JqeIRI0yW+geCjFy0G3ct2IgP1bKpWaicVfuifz94sc4vGnroOOJiBRZmSz0/8yfzM3DOzKh6kbScpN5qu4NtO/UR7NxEYlqhTrzxczONbPlZrbCzO4tYPt1ZpZhZvPDX10jH/XAZe3M5PFHz6HJ6DZ8lbyRF5I7MLPvJi67/lmVuYhEvX3O0M0sHngJOAtYA8wxs7HuviTf0Lfc/dZiyBgRX44byA2TurOoym4u3VGH/je+T52/pAUdS0QkYgqz5HICsMLdVwKY2SigPZC/0EulTeu/p9dz5zOw4lLqxccztuE/ueCaR4KOJSIScYUp9DrA6j2erwFOLGDcpWbWAvgW6OHuq/MPMLNuQDeA+vXrFz1tEXgoxNtDetB9xYtkVAjRY/dx9On1MSkH1SrW/YqIBKUwa+hWwGue7/mHQAN3PwaYBIwo6I3cfZC7p7l7Wo0aNYqWtAi+XziVtj0P5op1z1M3qzxzWr7Os4/NVZmLSEwrzAx9DVBvj+d1gXV7DnD3jXs8/Tfw5IFHK7rsXTt47pkO9N75CfHJ0L/CJdxy30jiyyUGEUdEpEQVptDnAEea2WHAWuAKoNOeA8ystruvDz+9EFga0ZSFMHP8YLpNuI1vquyi/Y5avNDtA+r99YSSjiEiEph9Frq755jZrcCnQDww1N0Xm1kfIN3dxwK3m9mFQA7wK3BdMWb+nS0bVtHrmbYMqLCYQxLiea/BvVx07eMltXsRkVLD3PMvh5eMtLQ0T09P3+/v91CIMcP+we3L+/FzhRC3ZTXlkbs+plK1QyKYUkSkdDGzue5e4DHXUXmm6A+LpnPrkEv5OHUDzbIrMPbsAaSd2TnoWCIigYq6Qv/g1fvotPxxLBmeLd+e254cRUJi+aBjiYgELupuetks7QLO31mXxVfNoMc976vMRUTCom6GXr/Rybz97B/OWRIRKfOiboYuIiIFU6GLiMQIFbqISIxQoYuIxAgVuohIjFChi4jECBW6iEiMUKGLiMQIFbqISIxQoYuIxAgVuohIjFChi4jECBW6iEiMUKGLiMQIFbqISIxQoYuIxAgVuohIjChUoZvZuWa23MxWmNm9BWxPMrO3wttnmVmDSAcVEZE/t89CN7N44CXgPKARcKWZNco37Hpgk7sfATwHPBnpoCIi8ucKc0/RE4AV7r4SwMxGAe2BJXuMaQ/0Dj8eDbxoZubuHsGsANxxB8yfH+l3FREpOU2bQr9+kX/fwiy51AH2vCvzmvBrBY5x9xxgC1At/xuZWTczSzez9IyMjP1LLCIiBSrMDN0KeC3/zLswY3D3QcAggLS0tP2avRfHv2oiIrGgMDP0NUC9PZ7XBdbtbYyZJQBVgF8jEVBERAqnMIU+BzjSzA4zs0TgCmBsvjFjgWvDjzsAk4tj/VxERPZun0su7p5jZrcCnwLxwFB3X2xmfYB0dx8LDAFeM7MV5M3MryjO0CIi8keFWUPH3ccB4/K99uAej3cBl0U2moiIFIXOFBURiREqdBGRGKFCFxGJESp0EZEYYUEdXWhmGcCP+/nt1YFfIhgnUkprLii92ZSraJSraGIx16HuXqOgDYEV+oEws3R3Tws6R36lNReU3mzKVTTKVTRlLZeWXEREYoQKXUQkRkRroQ8KOsBelNZcUHqzKVfRKFfRlKlcUbmGLiIifxStM3QREclHhS4iEiOiqtDNbKiZbTCzRUFn2ZOZ1TOzz81sqZktNrPuQWcCMLPyZjbbzBaEcz0cdKY9mVm8mX1tZh8FnT4akTQAAAORSURBVOU3ZvaDmX1jZvPNLD3oPL8xs1QzG21my8J/zk4uBZmOCv8+/fa11czuCDoXgJn1CP+ZX2RmI82sfNCZAMysezjT4uL4vYqqNXQzawFkAq+6e5Og8/zGzGoDtd19nplVAuYCF7n7kn18a3HnMqCiu2eaWTlgOtDd3WcGmes3ZnYnkAZUdvd2QeeBvEIH0ty9VJ2MYmYjgGnuPjh8X4Jkd98cdK7fhG8mvxY40d3394TBSGWpQ96f9UbuvtPM3gbGufvwgHM1AUaRd5/mLGA8cJO7fxepfUTVDN3dp1IK74Tk7uvdfV748TZgKX+872qJ8zyZ4aflwl+l4l9wM6sLnA8MDjpLaWdmlYEW5N13AHfPKk1lHtYG+E/QZb6HBKBC+A5qyfzxLmtBOBqY6e47wvdengJcHMkdRFWhRwMzawA0A2YFmyRPeFljPrABmOjupSIX0A+4GwgFHSQfByaY2Vwz6xZ0mLCGQAYwLLxENdjMKgYdKp8rgJFBhwBw97XA08AqYD2wxd0nBJsKgEVACzOrZmbJQFt+f3vPA6ZCjyAzSwHGAHe4+9ag8wC4e667NyXvXrAnhH/sC5SZtQM2uPvcoLMU4FR3Pw44D7glvMwXtATgOOAVd28GbAfuDTbS/4SXgC4E3gk6C4CZVQXaA4cBhwAVzezqYFOBuy8FngQmkrfcsgDIieQ+VOgREl6jHgO84e7vBp0nv/CP6F8A5wYcBeBU4MLwevUooLWZvR5spDzuvi786wbgPfLWO4O2Blizx09Xo8kr+NLiPGCeu/8cdJCwM4Hv3T3D3bOBd4FTAs4EgLsPcffj3L0FecvHEVs/BxV6RIQ/fBwCLHX3Z4PO8xszq2FmqeHHFcj7g74s2FTg7r3cva67NyDvR/XJ7h74DMrMKoY/1Ca8pHE2eT8mB8rdfwJWm9lR4ZfaAIF+4J7PlZSS5ZawVcBJZpYc/rvZhrzPtQJnZjXDv9YHLiHCv2+FuqdoaWFmI4FWQHUzWwM85O5Dgk0F5M04rwG+Ca9XA9wXvhdrkGoDI8JHIMQBb7t7qTlEsBQ6GHgvrwNIAN509/HBRvqv24A3wssbK4EuAecBILwWfBZwQ9BZfuPus8xsNDCPvCWNryk9lwAYY2bVgGzgFnffFMk3j6rDFkVEZO+05CIiEiNU6CIiMUKFLiISI1ToIiIxQoUuIhIjVOgiIjFChS4iEiP+H34gWdI3gR62AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "iter = 100\n",
    "range_x1 = np.linspace(1,9,17)\n",
    "results = []\n",
    "for sigma in range_x1:\n",
    "    calu = np.zeros(3)\n",
    "    for i in range(iter):\n",
    "        gussian_noise = np.random.normal(0,sigma,size=(num,1))\n",
    "        dectected_slope = slope_distance + gussian_noise\n",
    "        sol4_root = root(f5,[x_axis,y_axis,z_axis])#,method='hybr')\n",
    "        sol4_fsolve = fsolve(f5,[x_axis,y_axis,z_axis])\n",
    "        calu += (sol4_fsolve-gt)**2\n",
    "    calu /= iter\n",
    "    results.append(np.sqrt(calu))\n",
    "results = np.array(results).reshape(-1,3)\n",
    "plt.figure(1)\n",
    "plt.title(\"%d\" %(num))\n",
    "plt.plot(range_x1,results[:,0],'r',range_x1,results[:,1],'b',range_x1,results[:,2],'g')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.37084577, 0.        , 0.37081952],\n",
       "       [0.57119398, 0.        , 0.57118943],\n",
       "       [0.69938651, 0.        , 0.69937278],\n",
       "       [0.83812614, 0.        , 0.83817531],\n",
       "       [0.92297803, 0.        , 0.92299833],\n",
       "       [1.16716537, 0.        , 1.16715084],\n",
       "       [1.56187693, 0.        , 1.5618744 ],\n",
       "       [1.63753286, 0.        , 1.63752307],\n",
       "       [1.94840809, 0.        , 1.948404  ],\n",
       "       [2.0011161 , 0.        , 2.001074  ],\n",
       "       [2.05537186, 0.        , 2.05535391],\n",
       "       [2.2244824 , 0.        , 2.22448162],\n",
       "       [2.48512145, 0.        , 2.48514157],\n",
       "       [2.38039182, 0.        , 2.38040223],\n",
       "       [2.60731285, 0.        , 2.60730693],\n",
       "       [2.75396487, 0.        , 2.75398378],\n",
       "       [3.32399209, 0.        , 3.32399482]])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x20b79f49908>,\n",
       " <matplotlib.lines.Line2D at 0x20b79f50948>,\n",
       " <matplotlib.lines.Line2D at 0x20b79f50bc8>]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hUZfrG8e+ThCSEAEGKIEVEXVdgFTTYBQQrolgQFRVlZbGLiKuiqwJ2bGAFpNpABQsqIiBKU0pAQKqyqFQlIi20lHl+f2TcH8YgCUxyksn9ua5czMx5M+cW4ebNO6eYuyMiIqVfTNABREQkMlToIiJRQoUuIhIlVOgiIlFChS4iEiXigtpxtWrVvH79+kHtXkSkVJo7d+6v7l49v22BFXr9+vVJS0sLavciIqWSmf20t21achERiRIqdBGRKKFCFxGJEip0EZEooUIXEYkSKnQRkSihQhcRiRIqdBGRYtSndytmjh9cJO+tQhcRKSZfjRvIQ3zB+JlvFMn7q9BFRIpBKCeb7hN7UGt7DP/u9k6R7COwU/9FRMqStwd3Z3bKdoZVvZ4KVWoUyT40QxcRKWI7t/7GvSteocnm8nS6aUCR7UczdBGRItb/+Y6sSs5h2D8eJia26GpXhS4iUoQ2/LiYx3Z8xoU7a9Lqkh5Fui8tuYiIFKGHXunAzjjo23Foke9LhS4iUkQWf/UhgxKXcFPmMRzV7Lwi358KXUSkiNz1bhcqZRoP3Tq6WPanNXQRkSLw2duPMj7lV56Ob0vVOkcWyz41QxcRibDszF30mPUwDbbFcevtbxbbfjVDFxGJsKEvd2Fx5d2MrnsXCRUqFdt+NUMXEYmgbRvX8cD6tzhtUyUuue7JYt23ZugiIhH0RP/L2JDkfNS6HxZTvHNmFbqISISsWvI1z4a+4qrt9Tnh7M7Fvn8tuYiIRMh9QzoC8Nj1bwWy/30WupklmtlsM1tgZovNrHc+YxLM7G0zW2Fms8ysflGEFREpqWZPGMablX7kzphTqNfw5EAyFGSGvhto5e7HAk2Ac83spDxjrgc2ufsRwHNA8X4SICISIA+FuHNcN2rsMO7t9m5gOfZZ6J4rI/y0XPjL8wxrB4wIPx4NtDYzi1hKEZES7L3h9zCjyjYertWRilUPCSxHgdbQzSzWzOYDG4CJ7j4rz5DawGoAd88GtgBVIxlURKQk2r19K3cv6UfjLQn88+aiuVdoQRWo0N09x92bAHWAE8yscZ4h+c3G887iMbOuZpZmZmnp6emFTysiUsK8+PxVrKyYzTMnPURcfGKgWQp1lIu7bwa+BM7Ns2kNUBfAzOKAysBv+Xz/IHdPdffU6tWr71dgEZGS4tfVy3l468ecu7kaZ3foGXScAh3lUt3MUsKPywNnAsvyDBsLXBt+3B6Y7O5/mqGLiESTPi91YFs8PH1ZsEstvyvIiUW1gBFmFkvuPwDvuPvHZtYHSHP3scAQ4HUzW0HuzPyKIkssIlICLJs9jpcTFtJ1V0MandIu6DhAAQrd3RcCTfN5/cE9Hu8CLotsNBGRkuvukf8kqTz0vumdoKP8j079FxEppMnvPcNHKb/weOw51KjfKOg4/6NT/0VECiEnK5Me0x7g0G2x3HHHqKDj/IFm6CIihfDagBuZn7KTkbVuIzE5Jeg4f6AZuohIAW3ftIH7V4/gxM0VuLxLv6Dj/Ilm6CIiBfRU/w6srxBidLNniv1a5wWhQhcRKYC136XRN3sKHXbU5ZQ2NwQdJ18l758YEZES6D+DriTH4Inr3gg6yl6p0EVE9mHe5DcZkbyCbqFmHHZM86Dj7JUKXUTkL/y4aDq3fHQjVXcZ990e3LXOC0Jr6CIi+Zg3+U2e+uQ+3k1ehVWEoTVvIOXgQ4OO9ZdU6CIiYR4KMeHdx3nq62f4vMomKiZA9+xUunV6mTpHNQs63j6p0EWkzMvatYO3h/Xg6e+GsyBlF4fEx/Bk3HnccNMAKteoF3S8AlOhi0iZtW3jOga/ejPPbfyY1ck5NLQEhh7UmY539SOhQqWg4xWaCl1EypyfVy7k+WE38kr2TDYnOi2yKvPKEXdw3hX/ISa29NZi6U0uIlJIy2aP45l37+S1xOVkxcKl22vz7zMe5oSzOwcdLSJU6CIS9WZ88gp9J/VhbMrPJCbA9ZmNuPOK/hzRtHXQ0SJKhS4iUSmUk82Hr9/PU/Nf5usqGVRNMB705tzadQDV6x0ddLwioUIXkaiTNnEE1356A0sq7+awuDheSGpP59tfokKVGkFHK1IqdBGJGjlZmfR94gIezJpAzdhYRh1yO5de+yRx8YlBRysWKnQRiQqrFn/FNa+ex9QqW+mQUZcB/55ClVqHBR2rWOlaLiJS6o0aeBvHvHEq8ypsZUS1fzHq6R/LXJlDAQrdzOqa2RdmttTMFptZt3zGtDSzLWY2P/z1YNHEFRH5f1s3rKZTjwZc+fOLNNyZzIIOX9DplkEl8uYTxaEgSy7ZQA93n2dmFYG5ZjbR3ZfkGTfN3dtGPqKIyJ/N+Phlrv7ydlYn59CLltzf99Mys1a+N/ssdHdfD6wPP95mZkuB2kDeQhcRKXJZu3bw8OPn8CjTqU8c004ayMnndQ06VolQqJ9LzKw+0BSYlc/mk81sgZl9amaN9vL9Xc0szczS0tPTCx1WREqWbRvX8XCf1rzarxPpq5YW+f5WfPM5p99bg4djptNp+xHM7/mTynwP5u4FG2iWDEwBHnX39/JsqwSE3D3DzNoA/d39yL96v9TUVE9LS9vP2CIStJULvqTdiHNZVHk3ADEhaLElhfZ1z+Hii+6l1uFNIrYvD4UY/mIXbvtlGOVCxsDDu9Ghy3MRe//SxMzmuntqftsKNEM3s3LAGODNvGUO4O5b3T0j/HgcUM7Mqh1AZhEpwSa/9wzNRrZibUImExo+zoJW73I/p7M+Zge3ZLxN7deb0vyOyvTvewmrl80+oH39tnYFl91Vj39uGkaz7Sks7PR1mS3zfdnnDN3MDBgB/Obud+xlTE3gF3d3MzsBGA0c6n/x5pqhi5Q+Hgrx4jMd6J4xhr9vS+DDTuM4vEmrP4xZMnMsoz/rx+gtX/Nt5V0AnLQ5mfY1WnLpBfdQv/FpBd7f5DFP02nmPWwoH+KRhPPocfcHxJaLj+h/U2nzVzP0ghT6acA04FsgFH75PqAegLsPMLNbgZvIPSJmJ3Cnu3/1V++rQhcpXXZv38rNvVIZmvw97TbX5PX/zKVi1UP+8nu+S/uMMeOfZfTGacxL2QnA8ZuTaF+tOe3b3LXXi2Pt3r6VBx47k6fLzeFv2+J56/yhHNfqqoj/N5VGB1ToRUWFLlJ6/LxyIZc8fypfV8nggdDp9HpwcqGvG75ywZeM+eQpRqdPYXbKdgCO3ZxI+yqncum53Tn6hPMBWDrrYzqO6sD8lJ3cuKMhz9w3haTKWsH9nQpdRPZb2qTXuOizzmyKDzG8QXcuu/7ZA37PVUu+5r2P+jJ6/efMqLINgIZbEjg9rgEj4peSnG0MbXQ/F1zz8AHvK9qo0EVkv7w54Ga6rHmFg3fF8sEFb9KkxeUR38fa79J4/8MnGbN2ElMrb+bsrdUYdtvn1GxwTMT3FQ1U6CJSKDlZmfTsdRpPxc+hxabKvHvn18VyDfHd27eWynt5Fqe/KnRdbVFE/mDzLz/R8YlmfJqSzs07GtPviVmUS0wqln2rzA9M2byCjYjka9nscZzY90gmVkxnYOWreOnJb4utzOXAaYYuIgCMe6s3Vy7qRUKcMbnZC5x+wa1BR5JCUqGLlHEeCtH38fPpmTWeJjvL88G/Pqdew5ODjiX7QYUuUobt2PIrXfqkMrLST1y+tS5De83TMd+lmNbQRcqo1UtncXqveoyq+BOPxZ7NyKd/VJmXcpqhi5RBMz55hUun3sqO8iHGHvEAba/uE3QkiQAVukgZsmb5HHoPvoZhScs5LLscX1zyIUefqBuNRQsVukgZ8Ovq5Twx4CpejJmLJ8ItmcfS6673y+SNlKOZCl0kim3buI7nXryap3d/wfZycM32w+l13fBCXcJWSg8VukgU2r19KwNeuo5HN35AepJz8c5aPHLhABqedGHQ0aQIqdBFokhOViavD7yZh34azqrkHM7YncLjrZ7mxHOuDzqaFAMVukgU8FCID17ryf0L+7G0cibHZycxuOEDnHnp3ViMjk4uK1ToIqXc5Peeoee0h5idsp2/E8/oundxyQNPqsjLIBW6SCk1Z+II7vvkTiZV+Y26cbEMqXIdnXq+Qlx8YtDRJCAqdJFSZumsj3ng7RsZU3kt1RKN5xIv4sYew0hMTgk6mgRMhS5SSvy0eAa9h3dmRNL3JJWHXrSke4/XqVS9TtDRpIRQoYuUMNs3bWDp3PEsXjaVResXsHjHTyyK+43VyTnEJ0K3rOPoeeMbxXIHISld9lnoZlYXeA2oCYSAQe7eP88YA/oDbYAdwHXuPi/ycUWix66MzSyf+xmLlk5h8dr5LMpYyeLYjfyQnI1b7piEGDia8jQP1aVx7FF0vOwhXdpW9qogM/RsoIe7zzOzisBcM5vo7kv2GHMecGT460TglfCvImVe1s7tfPfNRBYv/pJFa+axeOt/WWzpfF8xi1D4QJQ4g6MsgVSvxXUxR9Ko7vE0/kdrGhzTQh9ySoHts9DdfT2wPvx4m5ktBWoDexZ6O+A1z73j9EwzSzGzWuHvFSkTMnds4/v5n7NkyRSWrPmGJVv/y2LS+S55N1mxuWNigCNiytHIq9PBDqdR7eNo3PAMjmzamvjyyYHml9KvUGvoZlYfaArMyrOpNrB6j+drwq/9odDNrCvQFaBevXqFSypSQuzK2Mx38yayZOlUlqybz+JtK1nCr3xfMZOc8IzbgAYx5Tg65yAu8AY0PrgJjY5uwd+PP0dHo0iRKXChm1kyMAa4w9235t2cz7f4n15wHwQMAkhNTf3TdpGSZMeWX1n+zUSWLJvGkrULWJyxkiUxG/lv8v8vlcQAR1g5Gno1LvHDaHRwExoe3Zyjmp5F+UoHBZpfyp4CFbqZlSO3zN909/fyGbIGqLvH8zrAugOPJ1K8fvlhEd1eaMOcmPV/+HAyzuBIi+cYr8GVdjgND8kt7r81PYuECpWCDS0SVpCjXAwYAix192f3MmwscKuZjSL3w9AtWj+X0mbtd2m0Hngqq8tn0nZnXa61I2hYpykNG7XgiGNbaY1bSryCzNBPBa4BvjWz+eHX7gPqAbj7AGAcuYcsriD3sMXOkY8qUnRWLfmaVkOasyEhm89OfonT2t4cdCSRQivIUS7TyX+NfM8xDtwSqVAixemHBVNo9XprNsXnMLHFYF1qVkotnSkqZdr3cyfQalQbdsSF+PzM1zm+9dVBRxLZbyp0KbOWzvyIVu9dRE6MM7nNKI5t3iHoSCIHRIUuZdLCaaM58+MOxGJ8efEHujWbRAUVupQ58ya/yVkTrqF8yJh8+Tj+lnpO0JFEIkK3NJEyZdZnQ2g98WoqZscw9erPVeYSVVToUmZM//hlzprShYOy4pjSeQoNjm0ZdCSRiFKhS5nwxfvPcu7Xt1Brdzmmdv2aQxudGnQkkYhToUvUm/D2Y7SZ24NDdyYw5ZY0av8tNehIIkVChS5R7eM3HuSCRfdz1PZEvrzjG2o2OCboSCJFRke5SNR6f9g9XP5DX47JSGLCvxdy0CGHBx1JpEip0CUqvT2oG1etfZ5m25IZ33MxlWvo+vsS/VToEnVee6krnTe8yqlbKvHJA0upWPWQoCOJFAsVukSVwf2vpeum1zhjSxXGPrSMClVqBB1JpNio0CVqvPR0B27d/i7nbqnGe32W645BUuao0KXUy9yZQf/nLufurHFcsPlg3n3kO91FSMokFbqUOlvT1/D15BFMXzKeaVu+ZVbyFnaVg0u31Oatx5bpzkJSZqnQpcT7eeVCpn8xgmkrPmf6ru+YX2knoZjckyiaWhI35jSl+RFncUHH3sTFJwYdVyQwKnQpUTwUYsU3nzNtxltMWzWd6f4TKyplAVA+Fk7yFO63VE5v2IaTWnXSESwie1ChS6CyM3exYNpopqWNYfovc5hebj2/JIUAOCjOOG33wdwQn8rpx11E0xZXEJ9YIeDEIiWXCl0C8fPKhfzrhbP4MmkDGfG5r9WPieOsnEM5vfKpnHZSB/5+QhtiYmKDDSpSiqjQpdjt3r6VS58/jW+St9E5uzGn1WrJ6S06UeeoZkFHEynV9lnoZjYUaAtscPfG+WxvCXwI/BB+6T137xPJkBI9PBTitt4n8lWVbbxd+w46dHku6EgiUaMgM/ThwIvAa38xZpq7t41IIolqA567ilcrLKNnzikqc5EI2+flc919KvBbMWSRKDd17AvcvmUU52+uwcMPfhl0HJGoE6nroZ9sZgvM7FMza7S3QWbW1czSzCwtPT09QruW0uCnxTNoP6Mbh2+P582ec4iNKxd0JJGoE4lCnwcc6u7HAi8AH+xtoLsPcvdUd0+tXr16BHYtpcGOLb9y8eCz2B3rfHjFB7qUrUgROeBCd/et7p4RfjwOKGdm1Q44mUQFD4W4vs/xzK+8k5GNe3FUs/OCjiQStQ640M2spplZ+PEJ4ffceKDvK9Gh7+PnM6rSKh6LO4c2HR8KOo5IVCvIYYsjgZZANTNbAzwElANw9wFAe+AmM8sGdgJXuLsXWWIpNca91ZueWeO5fGtd7nl6XNBxRKLePgvd3a/cx/YXyT2sUeR/vkv7jI7f9uLYXeUZ8tBcLEb3IxcpavpbJhG3ZcMq2r11IeXc+KDLRCqk6ANwkeKgU/8lokI52Vz9eDNWVMxk0vH9OLTRqUFHEikzVOgSUQ/2bsnHKRt4MbkDLdp1CzqOSJmiJReJmHeH3MmjsTPoknEUN3cfGXQckTJHM3SJiAVT3+G6lc9x8vZkXnx0tj4EFQmACl0O2K+rl9NubEeqWAxjbpumGzSLBESFLgcka9cOOjxzEj9XzGFa8+HUOrxJ0JFEyiwVuhyQHr1O5osqmxlR7V80O+vaoOOIlGla6JT9NvT5zrxQfiHddx9Hp1sGBR1HpMzTDF32y8zxg7kpfThnbjuIvn1nBB1HRFChy35Y9/08Lpl8A3VCcbx992zi4hODjiQiqNClkHZlbObil5uzNSnEhDajOeiQw4OOJCJhKnQpMA+FuKlXM2anbOe9Q++m8akXBx1JRPagD0WlQLZtXEfHfx/G8IoreNCbc/F1TwYdSUTy0Axd9mnRjPdp/94VfJ+cyaOxZ3Nvz0+CjiQi+VChy196/eUbuWHdQCrFxTCpybOccXH3oCOJyF6o0CVfuzI2c3ufk3m1wjKab6/MqFu/1FmgIiWcCl3+5L/zJ9N+RFvmp+zknuyTeOTJL3RookgpoEKXP3h/+D10Xt4XSzQ+OvwB2l7dJ+hIIlJAKnQBci+ydW+f5jybMJfUXUm888/xHPaP04OOJSKFsM/DFs1sqJltMLNFe9luZva8ma0ws4VmdlzkY0pRWvtdGmfcW5NnE+Zy847GTH9kvcpcpBQqyHHow4Fz/2L7ecCR4a+uwCsHHkuKy8R3n6DJ0BOYX2Ebb9W6lZee/FbXMxcppfa55OLuU82s/l8MaQe85u4OzDSzFDOr5e7rI5RRikBOViaPPHYOvf1Ljs5MYPTlozn6xLZBxxKRAxCJNfTawOo9nq8Jv/anQjezruTO4qlXr14Edi37I33VUq5+9nQmVNnI1RmHMeCBmVSoUiPoWCJygCJx6r/l85rnN9DdB7l7qrunVq9ePQK7lsL6atxAmr7QmCkVNzKw8tW81neFylwkSkRihr4GqLvH8zrAugi8r0SQh0L063sJd+/4kLoex1dnvMZxra4KOpaIRFAkZuhjgU7ho11OArZo/bxk2bJhFe3vqsuduz/k/IyazPv3CpW5SBTa5wzdzEYCLYFqZrYGeAgoB+DuA4BxQBtgBbAD6FxUYaXgsjN3kTb5DSbNHsWwrVP4qWI2T8W3pcfTH2IxusimSDQqyFEuV+5juwO3RCyR7BcPhVg2ZxyTpg7n83Uz+KL8z2xNyN12fE4SI5r157S2NwcbUkSKlM4ULcXWfT+PzycOYtLKiUyK+ZF1FUIAHBYbx+XZf+fMBudyxtldqV7v6ICTikhxUKGXIlvT1zBl/EAmLfmISbuXsaTybgCqxhmtM2tzZuWWtD7jehoc2zLYoCISCBV6CZa5M4OZE4fx+bzRTNr8DbMqbyMnBsrHwumhqlxXrjVnnnI1x55+GTGx+l8pUtapBUqY7KzdfPjGfxi6cARTktLZHg8xDs2owL1+Cmce056Tz/onCcmVg44qIiWMCr2E2LT+BwYPu40XN41nVXIO9eJiuS7nH5xZ/3xannsjKQcfGnREESnhVOgBWzZ7HM+PuZsRcYvZEQ8tsirTr94NXHhVH2LjE4KOJyKliAo9AKGcbD575zH6z36ez1I2khAPHXcewe1nP0KTFpcHHU9ESikVejHK+O1nRgzrxgvrPmB5pUxqlovhYWtF1+ufp0b9RkHHE5FSToVeDH5YOJUX376TIaG5bEmEZqEKvHHw9Vx2bV/iyycHHU9EooQKvYh4KMTUj16g35dPMrbSeiwO2m+rS7fm/+Gkc7ro9HsRiTgVeoTtytjMyGF30n/lSBak7OKgBOOe0MncfHV/6hzVLOh4IhLFVOgRsnrpLAaNuouBu2aQnuQ0sgQGpVzDVd2fI6lS1aDjiUgZoEI/ADlZmXz69iMMnDuQcZU24AZtMw+m2wl30eriO7WsIiLFSoW+H9Ysn8OQUXczOGMaa5JzqFkuhnv9FLpc+jiHHdM86HgiUkap0AsoJyuT8e88xsC0AXxS6RdCMXB2dlX61/snF1zZi3KJSUFHFJEyToW+D2uWz2HoqHsYnDGV1ck5HBwfwz1+Cl0ufowGx7QIOp6IyP+o0PORk5XJZ+8+zsA5A/i40s+EYuCs7Ko8V7czF3bsrdm4iJRIKvQ9rP0uLXc2vm0Kq5JzqBEfw92hk/nXJY/pGuMiUuKV+ULPycpkwugnGDj7FT6u9DM5MXBW1kE8U+c6LuzYW2dyikipUWYL3UMhRg66jftWDOSnijnUiDfuyjmRf138GIc3aRV0PBGRQiuThf7f+ZO5eXgHJlTZSGpOEk/VuYF2HftoNi4ipVqBznwxs3PNbLmZrTCze/PZfp2ZpZvZ/PBXl8hHPXCZOzN4/JFzaDy6NV8nbeSFpPbM7LuJy65/VmUuIqXePmfoZhYLvAScBawB5pjZWHdfkmfo2+5+axFkjIivxg3khkndWFR5N5fuqE3/Gz+g9t9Sg44lIhIxBVlyOQFY4e4rAcxsFNAOyFvoJdKm9T/Q87nzGVhhKXVjYxnb4D9ccM3DQccSEYm4ghR6bWD1Hs/XACfmM+5SM2sOfAd0d/fVeQeYWVegK0C9evUKn7YQPBTinSHd6bbiRdLLh+i++zj69PyE5INqFul+RUSCUpA1dMvnNc/z/COgvrsfA0wCRuT3Ru4+yN1T3T21evXqhUtaCD8snEqbHgdzxbrnqZOZyJwWb/DsY3NV5iIS1QoyQ18D1N3jeR1g3Z4D3H3jHk9fBZ488GiFl7VrB889055eOz8lNgn6l7+EW+4bSWy5+CDiiIgUq4IU+hzgSDM7DFgLXAF03HOAmdVy9/XhpxcCSyOasgBmjh9M1wm38W3lXbTbUZMXun5I3b+fUNwxREQCs89Cd/dsM7sV+AyIBYa6+2Iz6wOkuftY4HYzuxDIBn4DrivCzH+wZcMqej7ThgHlF3NIXCzv17+Xi659vLh2LyJSYph73uXw4pGamuppaWn7/f0eCjFm2L+5fXk/fikf4rbMJjx81ydUrHpIBFOKiJQsZjbX3fM95rpUnin646Lp3DrkUj5J2UDTrPKMPXsAqWd2CjqWiEigSl2hf/jafXRc/jiWBM8mtuO2J0cRF58YdCwRkcCVupteNk29gPN31mHxVTPofs8HKnMRkbBSN0Ov1/Bk3nn2T+csiYiUeaVuhi4iIvlToYuIRAkVuohIlFChi4hECRW6iEiUUKGLiEQJFbqISJRQoYuIRAkVuohIlFChi4hECRW6iEiUUKGLiEQJFbqISJRQoYuIRAkVuohIlFChi4hECRW6iEiUKFChm9m5ZrbczFaY2b35bE8ws7fD22eZWf1IBxURkb+2z0I3s1jgJeA8oCFwpZk1zDPsemCTux8BPAc8GemgIiLy1wpyT9ETgBXuvhLAzEYB7YAle4xpB/QKPx4NvGhm5u4ewawA3HEHzJ8f6XcVESk+TZpAv36Rf9+CLLnUBva8K/Oa8Gv5jnH3bGALUDXvG5lZVzNLM7O09PT0/UssIiL5KsgM3fJ5Le/MuyBjcPdBwCCA1NTU/Zq9F8W/aiIi0aAgM/Q1QN09ntcB1u1tjJnFAZWB3yIRUERECqYghT4HONLMDjOzeOAKYGyeMWOBa8OP2wOTi2L9XERE9m6fSy7unm1mtwKfAbHAUHdfbGZ9gDR3HwsMAV43sxXkzsyvKMrQIiLyZwVZQ8fdxwHj8rz24B6PdwGXRTaaiIgUhs4UFRGJEip0EZEooUIXEYkSKnQRkShhQR1daGbpwE/7+e3VgF8jGCdSSmouKLnZlKtwlKtwojHXoe5ePb8NgRX6gTCzNHdPDTpHXiU1F5TcbMpVOMpVOGUtl5ZcRESihApdRCRKlNZCHxR0gL0oqbmg5GZTrsJRrsIpU7lK5Rq6iIj8WWmdoYuISB4qdBGRKFGqCt3MhprZBjNbFHSWPZlZXTP7wsyWmtliM+sWdCYAM0s0s9lmtiCcq3fQmfZkZrFm9o2ZfRx0lt+Z2Y9m9q2ZzTeztKDz/M7MUsxstJktC/85O7kEZDoq/Pv0+9dWM7sj6FwAZtY9/Gd+kZmNNLPEoDMBmFm3cKbFRfF7VarW0M2sOZABvObujYPO8zszqwXUcvd5ZlYRmAtc5O5L9vGtRZ3LgArunmFm5YDpQDd3nxlkrt+Z2Z1AKlDJ3dsGnQdyCx1IdfcSdTKKmY0Aprn74PB9CZLcfXPQuX4Xvpn8WuBEd9/fEwYjlaU2uX/WG7r7TjN7Bxjn7sMDztUYGEXufV/dIVIAAALISURBVJozgfHATe7+faT2Uapm6O4+lRJ4JyR3X+/u88KPtwFL+fN9V4ud58oIPy0X/ioR/4KbWR3gfGBw0FlKOjOrBDQn974DuHtmSSrzsNbAf4Mu8z3EAeXDd1BL4s93WQvC0cBMd98RvvfyFODiSO6gVBV6aWBm9YGmwKxgk+QKL2vMBzYAE929ROQC+gF3A6Ggg+ThwAQzm2tmXYMOE9YASAeGhZeoBptZhaBD5XEFMDLoEADuvhZ4GlgFrAe2uPuEYFMBsAhobmZVzSwJaMMfb+95wFToEWRmycAY4A533xp0HgB3z3H3JuTeC/aE8I99gTKztsAGd58bdJZ8nOruxwHnAbeEl/mCFgccB7zi7k2B7cC9wUb6f+EloAuBd4POAmBmVYB2wGHAIUAFM7s62FTg7kuBJ4GJ5C63LACyI7kPFXqEhNeoxwBvuvt7QefJK/wj+pfAuQFHATgVuDC8Xj0KaGVmbwQbKZe7rwv/ugF4n9z1zqCtAdbs8dPVaHILvqQ4D5jn7r8EHSTsTOAHd0939yzgPeCUgDMB4O5D3P04d29O7vJxxNbPQYUeEeEPH4cAS9392aDz/M7MqptZSvhxeXL/oC8LNhW4e093r+Pu9cn9UX2yuwc+gzKzCuEPtQkvaZxN7o/JgXL3n4HVZnZU+KXWQKAfuOdxJSVkuSVsFXCSmSWF/262JvdzrcCZWY3wr/WAS4jw71uB7ilaUpjZSKAlUM3M1gAPufuQYFMBuTPOa4Bvw+vVAPeF78UapFrAiPARCDHAO+5eYg4RLIEOBt7P7QDigLfcfXywkf7nNuDN8PLGSqBzwHkACK8FnwXcEHSW37n7LDMbDcwjd0njG0rOJQDGmFlVIAu4xd03RfLNS9VhiyIisndachERiRIqdBGRKKFCFxGJEip0EZEooUIXEYkSKnQRkSihQhcRiRL/B50wM2QsN0gBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "iter = 100\n",
    "range_x1 = np.linspace(1,9,17)\n",
    "results = []\n",
    "for sigma in range_x1:\n",
    "    calu = np.zeros(3)\n",
    "    for i in range(iter):\n",
    "        gussian_noise = np.random.normal(0,sigma,size=(num,1))\n",
    "        dectected_slope = slope_distance + gussian_noise\n",
    "        sol4_root = root(f5,[x_axis,y_axis,z_axis])#,method='hybr')\n",
    "        sol4_fsolve = fsolve(f5,[x_axis,y_axis,z_axis])\n",
    "        calu += (sol4_fsolve-gt)**2\n",
    "    calu /= iter\n",
    "    results.append(np.sqrt(calu))\n",
    "results = np.array(results).reshape(-1,3)\n",
    "plt.figure(1)\n",
    "plt.plot(range_x1,results[:,0],'r',range_x1,results[:,1],'b',range_x1,results[:,2],'g')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.29770537, 0.        , 0.29769793],\n",
       "       [0.48703619, 0.        , 0.48701882],\n",
       "       [0.56746498, 0.        , 0.56749738],\n",
       "       [0.74591266, 0.        , 0.7459098 ],\n",
       "       [0.8295017 , 0.        , 0.82950123],\n",
       "       [1.12314019, 0.        , 1.12313439],\n",
       "       [1.1011735 , 0.        , 1.10113306],\n",
       "       [1.28848294, 0.        , 1.28847998],\n",
       "       [1.40248342, 0.        , 1.40251827],\n",
       "       [1.57523363, 0.        , 1.57522752],\n",
       "       [1.58804903, 0.        , 1.58805666],\n",
       "       [2.06828277, 0.        , 2.06830618],\n",
       "       [2.02486808, 0.        , 2.02487966],\n",
       "       [2.1310495 , 0.        , 2.13105431],\n",
       "       [2.17160721, 0.        , 2.17159802],\n",
       "       [2.38031159, 0.        , 2.38032129],\n",
       "       [2.6754533 , 0.        , 2.67547813]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results ####4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x20b7a280fc8>,\n",
       " <matplotlib.lines.Line2D at 0x20b7a3003c8>,\n",
       " <matplotlib.lines.Line2D at 0x20b7a300708>]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhU9dnG8e+TPTEbQlqR1QWtO9ooWvcFFaWuKFgVVJRqRVBwA0UWUVDcUFGLCEJRUVmUKipUcK9LQKkirlVkq4Q9gYQs87x/ZHxLMZAAk5yZyf25rrmYmXNyzp2Q3Pnld87MMXdHRERiX0LQAUREJDJU6CIicUKFLiISJ1ToIiJxQoUuIhInVOgiInFChS4iEidU6CKbMbO3zKzUzIrDt6+DziRSWyp0kV/r6e6Z4du+QYcRqS0VuohInFChi/zaMDNbaWbvm9kJQYcRqS3Te7mI/JeZtQO+BMqALsCjQFt3/z7QYCK1oEIX2QYzex141d0fCTqLSE005SKybQ5Y0CFEakOFLhJmZrlmdpqZpZlZkpldDBwHvBF0NpHaSAo6gEgUSQaGAr8DKoGvgHPcXeeiS0zQHLqISJzQlIuISJxQoYuIxAkVuohInFChi4jEicDOcmnSpIm3bt06qN2LiMSkuXPnrnT3vOqWBVborVu3pqCgIKjdi4jEJDNbtLVlmnIREYkTKnQRkTihQhcRiRMqdBGROKFCFxGJEyp0EZE4oUIXEYkTevtcEYk7ocoKxo26khXrl5OVmk1mek7VLXNXsjJ3JTOrMZk5eWTl/IbMRr8lPWtXLCH2x7cqdBGJK5s2rKfrgAN5IWdx1RMbw7dtSAhBZjlkViSQWZlIZiiZLJLJJJXMhDSyEtNplJJDXkYeedm70aRRM/LyWtHkt3uQ12wfsnZtGhW/EFToIhI31v68iHOGHcLbjdYxIuVMel43keI1P1O8rpCitT9TXLSKoqJVFG9YQ/HGtRSVrKV4UxHFZcUUeTHFvpFiL6GIUoq9jOVWRBFrKKKS1RairBxYFb5999/9plRAk00J5JWn0IQM8hIyaZKcS156Y/Iyf0uT3KbkNW5Jk7xW5O3ehsbN9iYpJS3in78KXUTiwuKFH9HhyeP5JnsTz/z2L/zp6lEApGXm0qTFvju9fQ+FKFq1jJXLvqPwP9+zcuViCtcsYWXRfygsL2Slr6GQ9RSygUUsp9B+Yh1Acfi25L/b6rPp99x/d+Tf+kSFLiIx7/P3ptDh5QspSgvxetv7OOm8vhHfhyUkkJ3XnOy85ux5yAm1+piykmJWLfuOlcu/p3DFj6xcvYTCtUs55Pe1+/jtpUIXkZg2Z9oDnPNJXzItgXc7vMjBx3YKOtL/S0nPpOlebWm6V9t62Z8KXURi1qTRvei2+BH2Lk3ltSvn0HL/o4KOFCgVuojEpAeGn03fTdM5tiibl2/5jEZN9wg6UuBU6CISU0KVFfS9/QgeSvuUTuua8behX5CWmRt0rKgQ/ImTIiK1VFq8li437cFDaZ/Su7Qtz4/4UWW+GRW6iMSENct/4LTbW/NizhLuS+nIg3fNJSFRkwyb01dDRKLe4oUfcfqTx/Nt9iaea3odXXo8HHSkqKRCF5Go9q93J9Ph750pTgvxxqH3c+K5fYKOFLVU6CIStX45xzzLE3jvjMkcdMz5QUeKajXOoZtZCzObY2YLzWyBmfWuZp0TzGydmX0Wvt1RN3FFpKF47q/XcdqnfWlRmso/r/hAZV4LtRmhVwB93X2emWUBc81slrt/ucV677p7x8hHFJGGxEMh7r/nLG4qe5Xj1ufw0i2f6hzzWqqx0N19ObA8fL/IzBYCzYAtC11EZKdUlpfR544jeDhtPhesa86EoZ/rtMTtsF2nLZpZa+BQ4KNqFh9lZvPN7DUzO2ArH9/DzArMrKCwsHC7w4pI/CpZv5out+zJw2nzub70UCaN+EFlvp1qXehmlglMAa539/VbLJ4HtHL3Q4BHgJeq24a7j3b3fHfPz8vL29HMIhJnln07j+MGtmBK9lLuTz2LB4fN0znmO6BWhW5myVSV+TPuPnXL5e6+3t2Lw/dnAMlm1iSiSUUkLs19cyKHP3k4CzM28tJe/elz68tBR4pZNf4KNDMDngIWuvsDW1lnN+Bnd3czO4KqXxSrIppUROLOi0/1odsPD5LniXxwenS99W0sqs3fNEcDlwKfm9ln4ef6Ay0B3P0JoBNwjZlVACVAF3f3OsgrInHAQyHuvPMUBjKHPxRnMbX3B/x2jwODjhXzanOWy3uA1bDOo8CjkQolIvGrZP1qrhh8KJOyf+LSoj0ZPXSuDn5GiN6cS0TqzbJv53H8wJY8n/UTw5NOZ/y936rMI0iHkUWkXsyb/QxnvdaVtRkhpu3Zj7O73h10pLijQheROjf5qb50/eEB8kjk/dOe55DjLgw6UlxSoYtInfFQiKFD23OHz+ao4kym9f6nDn7WIRW6iNQJHfysfzooKiIRt/z7z/7/4OewxNN08LOeaIQuIhG1+cHPqXvcyjndhgUdqcFQoYtIxEwZexOX/vs+mpDIe6dOou3xnYOO1KCo0EVkp3koxF1DT2WAv0m7Dbvw0nUfsNueBwcdq8FRoYvITilZv5rugw/juexFXLy+NWPu/FTz5QFRoYvIDite/R/a39mGD3OLuSuhPf1GvI4l6FyLoKjQRWSHeCjEZUPz+Ti7mBeb96FT9/uDjtTgqdBFZIcMu+t0puQsZUTKmSrzKKFCF5Ht9uozg7i9chYXFbWi74jpQceRMBW6iGyXbwre4E8LBnNISTpjBhVozjyK6H9CRGptfeESznn2LFJCxktXziIjR1eajCYaoYtIrYQqK+h69+F8k13GrLb30+qAo4OOJFtQoYtIrdw5tD0v5/6Hh9LP5cRz+wQdR6qhKRcRqdHLE/oziLfoWrQnvW6cHHQc2QqN0EVkmxZ+9AqXfjWM/JIMnhj8iQ6CRjH9z4jIVq39eRFnv3Ae6ZXG1B6zSc/eNehIsg0aoYtItSrLy7h4+OH8kFXO7MMfocV+7YKOJDVQoYtItQbeeRIzcgsZldmZY//YM+g4UguachGRX5ky9ibuSnyf7sX7cM0NzwYdR2pJI3QR+R9fvD+Nbt/fx5EbMxk1VAdBY4n+p0RixOKFH1FeurFO97F62fecPe1CssoTmPKXt0ndJbtO9yeRpUIXiQHP/fU6Wj1/JHsMyOLuO0+l8KeFEd9HZXkZF93bjsUZFUw98XF2b3NYxPchdavGQjezFmY2x8wWmtkCM+tdzTpmZg+b2Xdm9i8z03eCSIS89twQui59lCPXZrJ/eS63hWbR4sn9uezGvZn75sSI7af/4GOZ2WgVjzXuylEdekRsu1J/ajNCrwD6uvt+wJHAtWa2/xbrdADahG89gMcjmlKkgXr/1cc5f8FADipK57X+C5j50Cq+PP3vdC87kMmp35P/3qX84fosJo3utVPTMc8/2Zt7kz/m6o37c2Xv8RH8DKQ+1Vjo7r7c3eeF7xcBC4FmW6x2NjDBq3wI5JpZ04inFWlA/vXuZDq+dy0tSpJ5vfcn5PymJQD7tevIqHs+Z2nvRTyYdg4rEku5aPkjtLojiyGDT+LnH77Yrv3Mf+cFLl/0MEevyWLkoI/q4lORerJdc+hm1ho4FNjyf70ZsHizx0v4deljZj3MrMDMCgoLC7cvqUgD8v1nszn1lc5kViQw84o5/Kb1Ab9aJ+c3Lbn+lml8c28Jr7YZxCFluzKQObQcexCX9t2Tj2eOq3E/Kxd/zTnTL2bXTQlM7vUeKemZdfHpSD2pdaGbWSYwBbje3ddvubiaD/FfPeE+2t3z3T0/Ly9v+5KKNBDLvp1H+4mnUWHOzE4v1fg2tQmJSZzxp4G89lAhX50xgz+XH8zLqT/Q7p9X0O6GTJ554i+UlRT/6uMqykrpfP9RLE+vYNqpY9ltz4Pr6lOSelKrQjezZKrK/Bl3n1rNKkuAFps9bg4s2/l4Ig3L6mXfc9rjf6AwtYLXTh7Hfu06btfH73t4Bx4ePp8lNyzmkYxOrLMyLvn5cVoOymHgoONZ/v1n/7/uzYOOZnajNTzx2+4c3r5bpD8VCYC5/2og/b8rmBkwHljt7tdvZZ0zgZ7AGUA74GF3P2Jb283Pz/eCgoIdCi0SjzasWUH7wXsxN6uYGW3v5eTzb9rpbYYqK5g1+R4e+XAkM3IKSQxBpw0tOSCnDQP8TXqVHsLIYZ/VvCGJGmY2193zq11Wi0I/BngX+BwIhZ/uD7QEcPcnwqX/KHA6sBG43N232dYqdJH/Kisp5o/9WvOPnFW82Pomzrv83ojv47tP3+SxF27iKfuU9alw/JocZg1fRnJaRsT3JXVnpwq9rqjQRapUlpdx8S1783zOYsY06kb3Xk/X6f6KV/+HGVPvof2ZvWjUdI863ZdE3rYKXe/lIhIgD4XoefuhPJ+zmHuTz6jzMgfI3HU3LrzywTrfj9Q/vfRfJEADBh3HExlfckt5O27q/2rQcSTGaYQuEpAHh5/DXYnvc9WG3zFs+AdBx5E4oEIXCcD4UVfRZ9PLdFrXjMfvma+3qJWIUKGL1LOXJ/Sn+4oxtF+3KxPv+pLE5JSgI0mcUKGL1KO3XnqIzt8OI794F6besUDvNy4RpUIXqSdz35zIWR/fwF6lqbza91Myd90t6EgSZ1ToIvXgq49ncPqsrjSuSGLm1e/RuHmboCNJHNKRGJE69tOX/+TUF88iwWHmRTNotk+1rwkR2WkaoYvUocKfFnLqUyewLrWSt0+bRJvftw86ksQxFbpIHVm34ic6PJjPol3KmNnuEdoe3znoSBLnVOgidaBo1TI63H0A87M3Mm2fARz7x55BR5IGQIUuEmEb1qzgzDv34+OcYp5vdSMdLxkSdCRpIFToIhG0cd1K/jh4X97PWc+zzXtx/hUjgo4kDYgKXSRCSovXcs7AfXkrdy1/2+0aOl81MuhI0sDotEWRCNi0YT3nDdiHf+SuZmyT7lx89WNBR5IGSCN0kZ1UVlLMBbfvw2u5hTyZ25XLeo4JOpI0UBqhi+yE8tKNdOm/D3/P/ZnHMrtwZe/xQUeSBkyFLrKDKspKuaT/75iWu5yR6edxTd/ngo4kDZwKXWQHVJaX0a3f73ghZzH3pXSk181Tgo4kokIX2V6hygq699ufZ7MXMSzxNPr2+3vQkUQAFbrIdglVVvDnWw9kfNb3DOZEbr399aAjifw/FbpILXkoRM/+bRmT+TW3h47hjoGzg44k8j902qJILXgoxPW3HcbjGQu4pbwdQ4a8HXQkkV9RoYvUwEMhbhpwBA+nzeeGTYcxbOgHuqizRCUVusg2eChE/4FHc3/KXK4rOZj77/5EZS5RS4Uusg2DhpzI8KQPuXrj/owc9qnKXKKaCl1kK4YOOYUh9g7di/dh1PD5KnOJejV+h5rZWDNbYWZfbGX5CWa2zsw+C9/uiHxMkfp1z10dGOBv0rVoT0YPX0BCosY+Ev1q8136NPAoMGEb67zr7h0jkkhkO639eRH9HjiDNeVFpCQkk2JJpCQkk5qQQkpiMimJKZvdUklJSiE1OZ2UpFRSktP+e0tJJzU1g3fmT+cOn81F61sxdvhClbnEjBq/U939HTNrXfdRRHbMdfeeyHOZP7B3RQplhMI3Z5M5ZcAmA3egInzbVPM2L1jXnAnDvyIxOaVuw4tEUKSGHkeZ2XxgGXCjuy+obiUz6wH0AGjZsmWEdi0N2QtjbmBi9g8M4gQG3j9nq+tVlpdRVlLMppIiykqKKdu0kbLSDVX/btrIptINlJVtpKyshMTEZP5w+lUkpaTV42cisvPM3WteqWqE/oq7H1jNsmwg5O7FZnYGMNLd29S0zfz8fC8oKNj+xCJhS78p4KBxR9CmNIP371mpApYGwczmunt+dct2+rC9u6939+Lw/RlAspk12dntimxLqLKCyx87jU0JzsRLpqnMRYhAoZvZbmZm4ftHhLe5ame3K7Itox7owqxGq3mgyZ9o8/v2QccRiQo1zqGb2XPACUATM1sCDASSAdz9CaATcI2ZVQAlQBevzTyOyA5a+OHfuXndFM7c8Bt6DPhb0HFEokZtznK5qIblj1J1WqNInSsrKeaS5zuTmWqMuW6WXuwjshmdYCsxZciw05mXW8K01rey254HBx1HJKpoeCMx44MZf2WYvc8VxW04p9uwoOOIRB2N0CUmFK1axqWze9KKJB7q91bQcUSikgpdYkKf4Sfy4y4VvH34KLIa7x50HJGopCkXiXrTJ9zGmMxvuCV0FMd0/EvQcUSilkboEtV+/uELrlwwjEPL0hl098yg44hENY3QJWp5KMSVD5/M+hRnYudJpKRnBh1JJKpphC5Ra8zD3XgldwUPpZ/L/keeFXQckainEbpEpe8+fZMbCidyyppdua7vC0HHEYkJKnSJOhVlpVw64RySQ8a4a17XBSZEakk/KRJ1hg8/kw9zi5m0ey+a73t40HFEYoZG6BJVCmaNZ3DFbP60vhWdrxoZdByRmKIRukSNjetWcsnrPdgtIZFHb9z61YdEpHoqdIkaNw87ia+zy3jzoPto1HSPoOOIxBxNuUhUeH3SUEalf84Nmw7jpPP6Bh1HJCZphC6BW7XkW66YN5ADKlK5e8ibQccRiVkqdAmUh0Jc/cCJrMwM8dopE0jLzA06kkjMUqFLoCY+fg2Tc5ZyT1IHDjnuwqDjiMQ0FboEZtEX79Fz6WiO3ZhN3xEvBR1HJObpoKgEomTdKrqN6YgDE658lcTklKAjicQ8jdClXlVsKmHc4z0YvPRZljYKMb7JVbQ+8JigY4nEBRW61ItQZQVTxt3E7QtH8U12OUdVZPFM26Ecf3avoKOJxA0VutQpD4X4x+R76PfPoczN3cgBnsrLe97GHy8egiVoxk8kklToUmc+njmOfjP6MLvRWlolJTI+rwcX939E8+UidUSFLhG38KNXuP35PzM1Zxl5qcbIjPP5841jSd0lO+hoInFNhS4Rs/jLDxk0ritPZ3zLLmkw2E7khpsnktV496CjiTQIKnTZaSsXf82wxy9mVOJcPA16lx9Gv2ueIa/F74KOJtKg1FjoZjYW6AiscPcDq1luwEjgDGAjcJm7z4t0UIk+xav/w4OPXsyI0tlsSIZuG9sw6PLxtNz/qKCjiTRItRmhPw08CkzYyvIOQJvwrR3wePhfiVNlG4v466jLGbpyKisynHNLmjL0rCd0IWeRgNVY6O7+jpm13sYqZwMT3N2BD80s18yauvvyCGWUKFFeupFJT93AHT+M5cesCk7YlMvLJ93Hkad1DzqaiBCZOfRmwOLNHi8JP/erQjezHkAPgJYtW0Zg11If/j3/LZ6c3I9xZR/zc0aIQyvT+et+Q2jf6RadSy4SRSJR6FbNc17diu4+GhgNkJ+fX+06Eh3KSoqZ/uxARn/+NLMarSYhATqW7UaPg/9Mhy63k5Co4+ki0SYSP5VLgBabPW4OLIvAdiUA3336JmOm3Ma48o9ZkeG0TE5kiJ3E5X8aTvN9Dw86nohsQyQKfTrQ08wmUXUwdJ3mz2NLWUkxLz0zgNFfjOfNRmtITIQ/bmhKj7ZXc2qnW/XKTpEYUZvTFp8DTgCamNkSYCCQDODuTwAzqDpl8TuqTlu8vK7CSmR9O3cWT069jacrCijMcFolJXKnncwVl9zL7m0OCzqeiGyn2pzlclENyx24NmKJpE5t2rCel54dwF8XTGBOo7UkJsFZG5vS49C/0P78mzUaF4lhOrLVQHxT8AZPTrudpyvnsjLdaZ2UxF0J7bn8kntpulfboOOJSASo0OPc7Kn3c+c7Q3mr0VqSEuHskmb0OOxaTjn/Jp2pIhJn9BMdpwr+MYH+r9zArEaraZ6cyN2Jp3L5pSPYbc+Dg44mInVEhR5nvvrkNQY8dxWTc5bSOM14IPVsrun7NGmZuUFHE5E6pkKPE4sXfsTgsV0Zl/ENGekw0I+nT9+JZOc1DzqaiNQTFXqM2/Kta3uVHUr/a54hr+V+QUcTkXqmQo9RRauX8+CjF3Nf6Rw2JEPXDXsz6PKnaXXA0UFHE5GAqNBjzKYN66veunbVNAr11rUishkVeoyoLC9j4uieDPxhLIuyKjlxUy7DTrqPdnrrWhEJU6FHOQ+FePlvt3Hb/Af5MmcTv6/M4Mn9B3DK+TfrrWtF5H+o0KPYnGkP0u+dO/got5h9LJkXm/fh/AEjVOQiUi0VehT6/P2p3Di5BzNzV9E8KZExjbrRrd8TJKWkBR1NRKKYCj2KrF72PXc8fC6Pp35OTppxX0pH/jJgPOnZuwYdTURigAo9ClSWlzH6kW7cvuJ51qY6V5ceyJDrptK4eZugo4lIDFGhB+ztlx+m19u38K+cUk4ozWXkmU9y8LGdgo4lIjFIhR6QnxZ+yE1jLuSF7MW0TEzUAU8R2Wkq9HpWsn41Ix7qxPCyOXj4PVduvn0yGTlNgo4mIjFOhV5PPBRiyvhbuPHLB1mUWckFG5szovskvVRfRCJGhV4PPn9/Kr0nd2dO7loOqkxjziH3ccI51wcdS0TijAq9Dq1e/m8GjjyXx1L/RU6aMSqzMz36Pa3zyUWkTqjQ60BleRmjH72MAT9PYo1OQxSReqJCj7B3pj9Cr7duZn5OKceX5vDwmWN0GqKI1AsVegRUlpfx3owneOyd+3ghezEtEhN5ofkNdBpwn05DFJF6o0LfQWUlxcyePpKpcyfykn1NYYaTrtMQRSRAKvTtsHHdSt6YNoIpX7zAK8k/si4NMpPhzJIWnN+8Ex3OvZnMXXcLOqaINFAq9BqsW/ETr0y9h6nfvMRrGcsoSYZdk4zzKtpw3l5dOOWcPqRl5gYdU0REhV6dwp8W8vK04Uz98TX+kVVIeSI0TUrg8vIDOe/Qrhx35rUkp2UEHVNE5H/UqtDN7HRgJJAIjHH34VssvwwYASwNP/Wou4+JYM46t+TrT5g2/R6mLn2Td3LWEkqAPRKT6F2Zz3lHXUW7U68gIVG//0QketXYUGaWCIwC2gNLgE/MbLq7f7nFqs+7e886yFhnNqxZwWOPX87kwrf5OHcDAAckpHIbx3LeCddyyLEX6CwVEYkZtRlyHgF85+7/BjCzScDZwJaFHlO++/RNzv3bmXyRs4l8Mrg78VTOO+169j28Q9DRRER2SG0KvRmweLPHS4B21ax3vpkdB3wD3ODui6tZJyq8+swgLv5iMIkpxsz9h9H+gluDjiQistNqM59g1TznWzz+O9Da3Q8G/gGMr3ZDZj3MrMDMCgoLC7cvaQSEKisYPPhEOn43mD1L05l7yTsqcxGJG7Up9CVAi80eNweWbb6Cu69y903hh08Cv69uQ+4+2t3z3T0/Ly9vR/LusLU/L+KsG5sxiLfoVrQX7w9eQusDj6nXDCIidak2hf4J0MbM9jCzFKALMH3zFcys6WYPzwIWRi7izvv8vSnk37s3b2StYFRmZ8bd+40uvCwicafGOXR3rzCznsAbVJ22ONbdF5jZEKDA3acDvczsLKACWA1cVoeZt8uk0b3ovugRchITePvIJ/jDGX8OOpKISJ0w9y2nw+tHfn6+FxQU1Nn2K8pKuWXQMTyQOpej12Tx4nXv0HSvtnW2PxGR+mBmc909v7plcflKmRU/LqDzQ8fwVqO19Cw5iPuHfUBKembQsURE6lTcvWrm45nj+P2oQ/gwcy3jm1zFI8P/pTIXkQYhrkboY0Z249qVE9jdk/jglGc59ISLgo4kIlJv4qLQN21Yz3WD2/HkLl9xalFjnr3xn7rcm4g0ODFf6IsXfkSn0Sfzce4G+lf+gSEj5pCYnBJ0LBGRehfThf7WSw9x4T/7UJruTG11M+dedk/QkUREAhOThe6hEA/eey43l0ynTXkK086fxu+OOCPoWCIigYq5Qt+wZgXd78zn+ZzFnLu+KU/3/5jsvOZBxxIRCVzMnbb44jP9eSF7MXcnnsqU+5eozEVEwmJuhN7tL6Np+2572h7fOegoIiJRJeZG6JaQoDIXEalGzBW6iIhUT4UuIhInVOgiInFChS4iEidU6CIicUKFLiISJ1ToIiJxQoUuIhInVOgiInFChS4iEidU6CIicUKFLiISJ1ToIiJxQoUuIhInVOgiInFChS4iEidU6CIicaJWhW5mp5vZ12b2nZndWs3yVDN7Prz8IzNrHemgIiKybTUWupklAqOADsD+wEVmtv8Wq3UH1rj73sCDwD2RDioiIttWm4tEHwF85+7/BjCzScDZwJebrXM2MCh8fzLwqJmZu3sEswJw/fXw2WeR3qqISP1p2xYeeijy263NlEszYPFmj5eEn6t2HXevANYBjbfckJn1MLMCMysoLCzcscQiIlKt2ozQrZrnthx512Yd3H00MBogPz9/h0bvdfFbTUQkHtRmhL4EaLHZ4+bAsq2tY2ZJQA6wOhIBRUSkdmpT6J8AbcxsDzNLAboA07dYZzrQLXy/EzC7LubPRURk62qccnH3CjPrCbwBJAJj3X2BmQ0BCtx9OvAU8Dcz+46qkXmXugwtIiK/Vps5dNx9BjBji+fu2Ox+KXBBZKOJiMj20CtFRUTihApdRCROqNBFROKECl1EJE5YUGcXmlkhsGgHP7wJsDKCcSIlWnNB9GZTru2jXNsnHnO1cve86hYEVug7w8wK3D0/6BxbitZcEL3ZlGv7KNf2aWi5NOUiIhInVOgiInEiVgt9dNABtiJac0H0ZlOu7aNc26dB5YrJOXQREfm1WB2hi4jIFlToIiJxIqYK3czGmtkKM/si6CybM7MWZjbHzBaa2QIz6x10JgAzSzOzj81sfjjX4KAzbc7MEs3sUzN7JegsvzCzH83sczP7zMwKgs7zCzPLNbPJZvZV+PvsqCjItG/46/TLbb2ZXR90LgAzuyH8Pf+FmT1nZmlBZwIws97hTAvq4msVU3PoZnYcUAxMcPcDg87zCzNrCjR193lmlgXMBc5x9y9r+NC6zmXALu5ebGbJwHtAb3f/MMhcvzCzPkA+kO3uHYPOA1WFDuS7e1S9GMXMxgPvuvuY8HUJMtx9bdC5fhG+mPxSoJ277+gLBiOVpRlV3+v7u3uJmb0AzHD3pytwfJwAAALRSURBVAPOdSAwiarrNJcBrwPXuPu3kdpHTI3Q3f0dovBKSO6+3N3nhe8XAQv59XVX651XKQ4/TA7fouI3uJk1B84ExgSdJdqZWTZwHFXXHcDdy6KpzMNOBr4Pusw3kwSkh6+glsGvr7IWhP2AD919Y/jay28D50ZyBzFV6LHAzFoDhwIfBZukSnha4zNgBTDL3aMiF/AQcDMQCjrIFhyYaWZzzaxH0GHC9gQKgXHhKaoxZrZL0KG20AV4LugQAO6+FLgP+AlYDqxz95nBpgLgC+A4M2tsZhnAGfzv5T13mgo9gswsE5gCXO/u64POA+Dule7elqprwR4R/rMvUGbWEVjh7nODzlKNo939MKADcG14mi9oScBhwOPufiiwAbg12Ej/FZ4COgt4MegsAGbWCDgb2APYHdjFzC4JNhW4+0LgHmAWVdMt84GKSO5DhR4h4TnqKcAz7j416DxbCv+J/hZwesBRAI4GzgrPV08CTjKzicFGquLuy8L/rgCmUTXfGbQlwJLN/rqaTFXBR4sOwDx3/znoIGGnAD+4e6G7lwNTgT8EnAkAd3/K3Q9z9+Oomj6O2Pw5qNAjInzw8Slgobs/EHSeX5hZnpnlhu+nU/WN/lWwqcDd+7l7c3dvTdWf6rPdPfARlJntEj6oTXhK41Sq/kwOlLv/B1hsZvuGnzoZCPSA+xYuIkqmW8J+Ao40s4zwz+bJVB3XCpyZ/Sb8b0vgPCL8davVNUWjhZk9B5wANDGzJcBAd38q2FRA1YjzUuDz8Hw1QP/wtViD1BQYHz4DIQF4wd2j5hTBKPRbYFpVB5AEPOvurwcb6f9dBzwTnt74N3B5wHkACM8Ftwf+HHSWX7j7R2Y2GZhH1ZTGp0TPWwBMMbPGQDlwrbuvieTGY+q0RRER2TpNuYiIxAkVuohInFChi4jECRW6iEicUKGLiMQJFbqISJxQoYuIxIn/A+HwIxZgAaZZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "iter = 100\n",
    "range_x1 = np.linspace(1,9,17)\n",
    "results = []\n",
    "for sigma in range_x1:\n",
    "    calu = np.zeros(3)\n",
    "    for i in range(iter):\n",
    "        gussian_noise = np.random.normal(0,sigma,size=(num,1))\n",
    "        dectected_slope = slope_distance + gussian_noise\n",
    "        sol4_root = root(f5,[x_axis,y_axis,z_axis])#,method='hybr')\n",
    "        sol4_fsolve = fsolve(f5,[x_axis,y_axis,z_axis])\n",
    "        calu += (sol4_fsolve-gt)**2\n",
    "    calu /= iter\n",
    "    results.append(np.sqrt(calu))\n",
    "results = np.array(results).reshape(-1,3)\n",
    "plt.figure(1)\n",
    "plt.title(\"%d\" %(num))\n",
    "plt.plot(range_x1,results[:,0],'r',range_x1,results[:,1],'b',range_x1,results[:,2],'g')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.30990996, 0.        , 0.30991397],\n",
       "       [0.51641654, 0.        , 0.51641679],\n",
       "       [0.60039239, 0.        , 0.60042056],\n",
       "       [0.75299696, 0.        , 0.75300848],\n",
       "       [0.8915473 , 0.        , 0.89153156],\n",
       "       [0.99579858, 0.        , 0.9957915 ],\n",
       "       [1.13481181, 0.        , 1.13478367],\n",
       "       [1.50994607, 0.        , 1.50992315],\n",
       "       [1.50365299, 0.        , 1.50364546],\n",
       "       [1.77070908, 0.        , 1.77070847],\n",
       "       [2.05323799, 0.        , 2.05327385],\n",
       "       [1.94164791, 0.        , 1.94165919],\n",
       "       [2.26691584, 0.        , 2.26691862],\n",
       "       [2.48499527, 0.        , 2.48498591],\n",
       "       [2.71726379, 0.        , 2.71727408],\n",
       "       [2.70465684, 0.        , 2.70467874],\n",
       "       [2.68369105, 0.        , 2.68367369]])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iter = 100\n",
    "range_x1 = np.linspace(1,9,17)\n",
    "results = []\n",
    "for sigma in range_x1:\n",
    "    calu = np.zeros(3)\n",
    "    for i in range(iter):\n",
    "        gussian_noise = np.random.normal(0,sigma,size=(num,1))\n",
    "        dectected_slope = slope_distance + gussian_noise\n",
    "        sol4_root = root(f5,[x_axis,y_axis,z_axis])#,method='hybr')\n",
    "        sol4_fsolve = fsolve(f5,[x_axis,y_axis,z_axis])\n",
    "        calu += (sol4_fsolve-gt)**2\n",
    "    calu /= iter\n",
    "    results.append(np.sqrt(calu))\n",
    "results = np.array(results).reshape(-1,3)\n",
    "plt.figure(1)\n",
    "plt.title(\"%d\" %(num))\n",
    "plt.plot(range_x1,results[:,0],'r',range_x1,results[:,1],'b',range_x1,results[:,2],'g')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f5(x):\n",
    "    m = np.array((x[0],x[1],x[2]),dtype = np.float64).reshape(1,3)\n",
    "    minus = (m - Paxis)\n",
    "    k = np.concatenate((minus**2, -distL1**2),axis = 1)\n",
    "    l = np.sum(k,keepdims = True,axis =1).reshape(-1,1)\n",
    "    return np.sum(l * minus,axis = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "Paxis = np.array((0,0,1,0,0,-1,1,0,0,-1,0,0)).reshape(4,3)\n",
    "distL1 = np.array((1,1,1,1)).reshape(4,1)\n",
    "sol4_root = root(f5,[233,323,233])#,method='hybr')\n",
    "sol4_fsolve = fsolve(f5,[233,323,233])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-3.63750972e-17,  0.00000000e+00, -3.47370894e-27])"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol4_fsolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
